Tugas Praktikum PSD 1 (Pertemuan 8 dan 9)
Dhea Puspita Adinda - G1401201090
2022-11-03
Nonlinear Regression Part 1
Library
library("MultiKink")
library("MatrixModels")
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(purrr)
library(rsample)Data Bangkitan
Data yang digunakan merupakan data bangkitan yang terdiri dari 1000 data, dengan mean sebesar 1 dan standar deviasi sebesar 1.
set.seed(123)
data.x <- rnorm(1000,1,1) #1000 data, miu 1, st dev 1
err <- rnorm(1000) #miu 0 st, dev 1
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))
Pada gambar di atas terlihat plot tersebut tidak membentuk garis lurus (tidak linear).
Regresi linier
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")
summary(lin.mod)##
## Call:
## lm(formula = y ~ data.x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.686 -2.574 -1.428 1.195 27.185
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6056 0.1902 24.22 <2e-16 ***
## data.x 8.3790 0.1340 62.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.2 on 998 degrees of freedom
## Multiple R-squared: 0.7967, Adjusted R-squared: 0.7965
## F-statistic: 3911 on 1 and 998 DF, p-value: < 2.2e-16Pada gambar di atas terlihat garis regresi linear merah tidak mendekati nilai observasi. Fungsi ini kurang tepat jika menggunakan fungsi linear biasa. Didapatkan Adjusted R-squared sebesar 0.7965 yang artinya sebanyak 79% keragaman data dapat dijelaskan oleh model.
Polynomial Ordo 2
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)
summary(pol.mod)##
## Call:
## lm(formula = y ~ data.x + I(data.x^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0319 -0.6942 0.0049 0.7116 3.2855
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.95193 0.04568 108.41 <2e-16 ***
## data.x 2.10732 0.05861 35.95 <2e-16 ***
## I(data.x^2) 2.99081 0.02338 127.93 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.007 on 997 degrees of freedom
## Multiple R-squared: 0.9883, Adjusted R-squared: 0.9883
## F-statistic: 4.221e+04 on 2 and 997 DF, p-value: < 2.2e-16Pada gambar di atas terlihat garis biru semakin mendekati nilai observasi, sehingga error semakin kecil. Fungsi ini lebih tepat jika menggunakan Polynomial dibanding menggunakan fungsi Regresi Linear. Didapatkan Adjusted R-squared sebesar 0.9883 yang artinya sebanyak 98% keragaman data dapat dijelaskan oleh model.
Fungsi Tangga
#regresi fungsi tangga
range(data.x)## [1] -1.809775 4.241040c1 <- 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))
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")
summary(step.mod)##
## Call:
## lm(formula = y ~ c1 + c2 + c3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.395 -3.534 -0.530 2.527 36.876
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.4499 0.4087 74.50 <2e-16 ***
## c11 -25.2395 0.5710 -44.21 <2e-16 ***
## c21 -19.4184 0.4536 -42.81 <2e-16 ***
## c31 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.121 on 997 degrees of freedom
## Multiple R-squared: 0.698, Adjusted R-squared: 0.6974
## F-statistic: 1152 on 2 and 997 DF, p-value: < 2.2e-16Pada gambar di atas, garis hijau membentuk seperti tangga. Garis tersebut terbagi menjadi 3 bagian yang mendekati nilai-nilai observasi pada bagian tersebut. Bagian 1 ada pada rentang -1.809775 sampai dengan 0, bagian 2 ada pada rentang 0 sampai dengan 2, bagian 3 ada pada rentang 2 sampai dengan 4.241040. Didapatkan Adjusted R-squared sebesar 0.6974 yang artinya sebanyak 69% keragaman data dapat dijelaskan oleh model.
Komparasi 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 5711.836
## 2 Poly (ordo=2) 2856.470
## 3 Tangga (breaks=3) 6109.609AIC yang terendah menandakan model yang terbaik, pada output di atas, didapatkan AIC terendah sebesar 2856.470 pada Model Polynomial Ordo = 2, maka model inilah yang terbaik digunakan.
MSE = function(pred,actual){
mean((pred-actual)^2)
}
MSE(predict(lin.mod),y)## [1] 17.60106MSE(predict(pol.mod),y)## [1] 1.010649MSE(predict(step.mod),y)## [1] 26.14693MSE yang terendah menandakan model yang terbaik, pada output di atas, didapatkan MSE terendah sebesar 1.010649 pada Model Polynomial Ordo = 2, maka model inilah yang terbaik digunakan.
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")Scatterplot
ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
theme_bw() 
Pada scatterplot di atas terlihat bahwa plot tidak linear. Akan dicari model yang bisa merepresentasikan pola hubungan tersebut.
Regresi Linear
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.502662583Pada output di atas, terlihat nilai Adjusted R-squared yang cukup kecil yaitu sebesar 0.3365 yang artinya sebanyak 33% keragaman data dapat dijelaskan oleh model.
ggplot(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()
Pada gambar di atas terlihat garis regresi linear biru tidak terlalu mendekati nilai observasi. Nilai observasi membentuk pola sedikit melengkung, sedangkan garis regresi linear merupakan garis lurus. Fungsi ini kurang tepat jika menggunakan fungsi Regresi Linear. Didapatkan Adjusted R-squared sebesar 0.7965 yang artinya sebanyak 79% keragaman data dapat dijelaskan oleh model.
Regresi Polinomial 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.89Berdasarkan output di atas, didapatkan Adjusted R-squared sebesar 0.3362 yang artinya sebanyak 33% keragaman data dapat dijelaskan oleh model. Didapatkan pula nilai AIC sebesar 5012.89.
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 = 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()
Pada gambar di atas, mirip seperti plot pada fungsi Regresi Linear, terlihat garis biru tidak terlalu mendekati nilai observasi.
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()
Berdasarkan output di atas, didapatkan Adjusted R-squared sebesar 0.3815 yang artinya sebanyak 38% keragaman data dapat dijelaskan oleh model. Pada gambar di atas, terlihat garis biru melengkung mendekati nilai observasi.
AIC <- cbind(AIC(mod_linear),AIC(mod_polinomial2),AIC(mod_polinomial3))
AIC## [,1] [,2] [,3]
## [1,] 5011.515 5012.89 4950.774MSE(predict(mod_linear),triceps$triceps)## [1] 16.01758MSE(predict(mod_polinomial2),triceps$triceps)## [1] 16.00636MSE(predict(mod_polinomial3),triceps$triceps)## [1] 14.89621AIC dan MSE yang terendah menandakan model yang terbaik, pada output di atas, didapatkan MSE terendah sebesar 14.89621 dan AIC terendah sebesar 4950.774 pada Model Polynomial Ordo = 3, maka model inilah yang terbaik digunakan.
Regresi 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-16Berdasarkan output di atas, didapatkan Adjusted R-squared sebesar 0.3588 yang artinya sebanyak 35% keragaman data dapat dijelaskan oleh model.
ggplot(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()
Pada gambar di atas, garis biru membentuk seperti tangga. Garis tersebut terbagi menjadi 5 bagian yang mendekati nilai-nilai observasi pada bagian tersebut.
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-16Berdasarkan output di atas, didapatkan Adjusted R-squared sebesar 0.4014 yang artinya sebanyak 40% keragaman data dapat dijelaskan oleh model.
ggplot(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()
Pada gambar di atas, garis biru membentuk seperti tangga. Garis tersebut terbagi menjadi 7 bagian yang mendekati nilai-nilai observasi pada bagian tersebut.
Perbandingan Model
MSE = function(pred,actual){
mean((pred-actual)^2)
}
# Membandingkan model
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.36779MSE yang terendah menandakan model yang terbaik, pada output di atas, didapatkan MSE terendah sebesar 14.36779 pada Model Tangga (breaks=7), maka model inilah yang terbaik digunakan.
Evaluasi Model Menggunakan Cross Validation
Regresi Linear
library(mlr3measures)## In order to avoid name clashes, do not attach 'mlr3measures'. Instead, only load the namespace with `requireNamespace("mlrmeasures")` and access the measures directly via `::`, e.g. `mlr3measures::auc()`.set.seed(123)
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 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.65 2.82
## 2 4.62 3.22
## 3 4.38 3.00
## 4 3.85 2.80
## 5 3.08 2.36
## 6 3.83 2.81
## 7 3.59 2.78
## 8 4.66 3.06
## 9 3.50 2.56
## 10 4.59 2.93# Menghitung Rata-Rata untuk 10 Folds
mean_metric_linear <- colMeans(metric_linear)
mean_metric_linear## rmse mae
## 3.973249 2.833886Polynomial Derajat 2 (ordo 2)
set.seed(123)
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 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.64 2.82
## 2 4.62 3.26
## 3 4.39 3.02
## 4 3.85 2.81
## 5 3.10 2.38
## 6 3.82 2.81
## 7 3.62 2.83
## 8 4.65 3.07
## 9 3.50 2.57
## 10 4.59 2.94# Menghitung Rata-Rata untuk 10 Folds
mean_metric_poly2 <- colMeans(metric_poly2)
mean_metric_poly2## rmse mae
## 3.977777 2.851787Polynomial Derajat 3 (ordo 3)
set.seed(123)
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 x 2
## rmse mae
## <dbl> <dbl>
## 1 3.49 2.60
## 2 4.48 2.99
## 3 4.21 2.85
## 4 4.02 2.75
## 5 3.03 2.09
## 6 3.63 2.59
## 7 3.53 2.52
## 8 4.54 2.91
## 9 3.27 2.33
## 10 4.27 2.68# Menghitung Rata-Rata untuk 10 Folds
mean_metric_poly3 <- colMeans(metric_poly3)
mean_metric_poly3## rmse mae
## 3.845976 2.632125Fungsi Tangga
set.seed(123)
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.835357 2.618775
## 2 4 3.882932 2.651911
## 3 5 3.917840 2.724368
## 4 6 3.836068 2.622939
## 5 7 3.789715 2.555062
## 6 8 3.812789 2.555563
## 7 9 3.781720 2.518706
## 8 10 3.795877 2.529479#berdasarkan rmse
best_tangga %>% slice_min(rmse)## breaks rmse mae
## 1 9 3.78172 2.518706#berdasarkan mae
best_tangga %>% slice_min(mae)## breaks rmse mae
## 1 9 3.78172 2.518706#Perbandingan 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.973249 2.833886
## 2 Poly2 3.977777 2.851787
## 3 Poly3 3.845976 2.632125
## 4 Tangga (breaks=9) 3.781720 2.518706RMSE dan MAE yang terendah menandakan model yang terbaik, pada output di atas, didapatkan RMSE terendah sebesar 3.781720 dan MAE terendah sebesar 2.518706 pada Model Tangga (breaks=9), maka model inilah yang terbaik digunakan.
Nonlinear Regression Part 2
Library
library(ggplot2)
library(dplyr)
library(purrr)
library(rsample)
library(splines)Data Bangkitan
Data yang digunakan merupakan data bangkitan yang terdiri dari 1000 data, dengan mean sebesar 1 dan standar deviasi sebesar 1.
set.seed(123)
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(-2,5), ylim=c(-10,70), pch=16, col="orange")
abline(v=1, col="red", lty=2)
Piecewise cubic polynomial (Fungsi Tangga)
dt.all <- cbind(y,data.x)
##knots=1
dt1 <- dt.all[data.x <=1,]
dim(dt1)## [1] 495 2dt2 <- dt.all[data.x >1,]
dim(dt2)## [1] 505 2Dari output di atas terlihat bahwa, dengan knot=1, tedapat 495 dimensi data pada data yang nilainya kurang dari 1, dan 505 dimensi data pada data yang nilainya lebih dari 1.
plot(data.x,y,xlim=c(-2,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)
Pada plot di atas terlihat bahwa terdapat garis biru yang terputus, yang merupakan 2 fungsi. Masing-masing fungsi dipisahkan oleh knot. Jarak yang putus/belum tersambung akan disambung dengan menggunakan Truncated Power Basis.
Truncated Power Basis
#dengan menggunakan truncated power basis
plot(data.x,y,xlim=c( -2,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 di atas terlihat bahwa terdapat garis hijau yang telah tersambung.
Perbandingan dengan k-fold CV
plot( data.x,y,xlim=c(-2,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)
Pada plot di atas terlihat bahwa terdapat garis hijau yang terputus, yang merupakan 3 fungsi. Masing-masing fungsi dipisahkan oleh knot=2 yang memisahkan ada di 0 dan 2. Jarak yang putus/belum tersambung akan disambung dengan menggunakan Truncated Power Basis.
Perbandingan 1 knot dan 2 knot dengan 5-fold CV
##1 knot
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] 22.08767 21.39598 29.47677 29.68683 25.18070mean(res)## [1] 25.56559Dari output di atas dapat dilihat bahwa rata-rata residu knot=1 adalah 25.56559
##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] 22.32898 21.33868 29.40459 29.79913 25.22047mean(res2)## [1] 25.61837Dari output di atas dapat dilihat bahwa rata-rata residu knot=2 adalah 25.61837
pada lipatan 1 knot = 1 lebih baik karena lebih kecil pada lipatan 2 knot = 2 lebih baik karena lebih kecil pada lipatan 3 knot = 2 lebih baik karena lebih kecil pada lipatan 4 knot = 1 lebih baik karena lebih kecil pada lipatan 5 knot = 1 lebih baik karena lebih kecil
Rata-rata residu pada knot = 1 lebih kecil daripada knot = 2, maka knot = 1 lebih baik.
Smoothing Spline
ss1 <- smooth.spline(data.x,y,all.knots = T)
ss1## Call:
## smooth.spline(x = data.x, y = y, all.knots = T)
##
## Smoothing Parameter spar= 1.499938 lambda= 0.0002745903 (26 iterations)
## Equivalent Degrees of Freedom (Df): 14.84845
## Penalized Criterion (RSS): 24848.46
## GCV: 211.964plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70), pch=16,col="orange")
lines(ss1,col="purple", lwd=2)
Berdasarkan output di atas, didapatkan spar = 1.499938 dan lambda optimum = 0.0002745903 dengan 26 iterasi, df = 14.84845, Penalized Criterion (RSS) = 24848.46, dan GCV = 211.964
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")Scatterplot
ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
theme_bw()
Berdasarkan pola hubungan yang terlihat pada scatterplot, kita akan mencoba untuk mencari model yang bisa merepresentasikan pola hubungan tersebut dengan baik.
Regresi Spline
#Menentukan banyaknya fungsi basis dan knots
dim(bs(triceps$age, knots = c(10, 20,40)))## [1] 892 6#dipilih 10,20,40 karena di rentang 0-10 terlihat plot membentuk suatu model, di rentang 10-20 terlihat plot membentuk suatu model, dst.
#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-16data_hasil <- cbind("y_aktual"=triceps$triceps,
"y_pred"=mod_spline_1$fitted.values,
"residual"=resid(mod_spline_1))
data_hasil## y_aktual y_pred residual
## 1 3.4 7.333196 -3.933196e+00
## 2 4.0 6.380136 -2.380136e+00
## 3 4.2 6.426062 -2.226062e+00
## 4 4.2 7.052817 -2.852817e+00
## 5 4.4 6.455335 -2.055335e+00
## 6 4.4 7.241677 -2.841677e+00
## 7 4.8 6.780652 -1.980652e+00
## 8 4.8 6.422456 -1.622456e+00
## 9 4.8 6.830448 -2.030448e+00
## 10 5.0 6.411711 -1.411711e+00
## 11 5.0 6.989901 -1.989901e+00
## 12 5.0 6.999516 -1.999516e+00
## 13 5.2 6.825882 -1.625882e+00
## 14 5.2 6.408154 -1.208154e+00
## 15 5.2 6.397555 -1.197555e+00
## 16 5.2 7.272060 -2.072061e+00
## 17 5.2 6.942197 -1.742197e+00
## 18 5.2 6.440604 -1.240605e+00
## 19 5.4 6.408154 -1.008154e+00
## 20 5.4 7.353674 -1.953674e+00
## 21 5.4 6.563489 -1.163489e+00
## 22 5.4 6.411711 -1.011711e+00
## 23 5.4 6.825882 -1.425882e+00
## 24 5.4 7.231580 -1.831580e+00
## 25 5.4 6.839602 -1.439602e+00
## 26 5.4 6.563489 -1.163489e+00
## 27 5.5 7.009155 -1.509155e+00
## 28 5.6 7.141440 -1.541440e+00
## 29 5.6 7.111718 -1.511718e+00
## 30 5.6 6.426062 -8.260622e-01
## 31 5.6 6.387066 -7.870665e-01
## 32 5.6 6.666944 -1.066945e+00
## 33 5.8 7.379334 -1.579334e+00
## 34 5.8 6.608254 -8.082537e-01
## 35 5.8 6.633177 -8.331771e-01
## 36 5.8 6.604134 -8.041339e-01
## 37 5.8 7.410211 -1.610211e+00
## 38 6.0 7.297476 -1.297476e+00
## 39 6.0 6.662691 -6.626910e-01
## 40 6.0 7.106781 -1.106781e+00
## 41 6.0 6.527793 -5.277934e-01
## 42 6.0 6.645770 -6.457696e-01
## 43 6.0 6.965970 -9.659704e-01
## 44 6.0 6.551496 -5.514959e-01
## 45 6.2 7.111718 -9.117177e-01
## 46 6.2 6.999516 -7.995164e-01
## 47 6.2 6.579625 -3.796256e-01
## 48 6.2 6.913884 -7.138837e-01
## 49 6.4 7.302570 -9.025696e-01
## 50 6.4 7.028503 -6.285031e-01
## 51 6.4 6.459047 -5.904694e-02
## 52 6.6 6.595924 4.075613e-03
## 53 6.6 7.028503 -4.285033e-01
## 54 6.6 7.009155 -4.091550e-01
## 55 6.6 7.038212 -4.382121e-01
## 56 6.8 7.251791 -4.517905e-01
## 57 6.8 6.980310 -1.803101e-01
## 58 6.8 7.410211 -6.102107e-01
## 59 7.0 6.666944 3.330556e-01
## 60 7.0 6.604134 3.958659e-01
## 61 7.0 6.612384 3.876165e-01
## 62 7.0 6.994706 5.294490e-03
## 63 7.0 6.895143 1.048568e-01
## 64 7.0 7.297476 -2.974763e-01
## 65 7.0 6.881161 1.188390e-01
## 66 7.0 6.684051 3.159489e-01
## 67 7.2 6.600024 5.999756e-01
## 68 7.2 7.176328 2.367131e-02
## 69 7.4 7.374197 2.580319e-02
## 70 7.4 6.512206 8.877939e-01
## 71 8.0 6.349577 1.650423e+00
## 72 8.0 7.307666 6.923337e-01
## 73 8.4 6.994706 1.405294e+00
## 74 8.6 7.067471 1.532529e+00
## 75 8.8 7.312766 1.487234e+00
## 76 8.8 6.380136 2.419864e+00
## 77 9.0 7.231580 1.768420e+00
## 78 9.0 7.405059 1.594941e+00
## 79 9.8 7.420522 2.379478e+00
## 80 9.8 7.092001 2.707999e+00
## 81 10.0 6.932732 3.067268e+00
## 82 10.2 6.975524 3.224476e+00
## 83 10.2 6.771702 3.428298e+00
## 84 13.0 6.624831 6.375169e+00
## 85 14.2 6.583685 7.616315e+00
## 86 16.0 6.951687 9.048313e+00
## 87 9.4 14.606604 -5.206604e+00
## 88 8.2 8.448493 -2.484932e-01
## 89 21.2 14.694642 6.505359e+00
## 90 9.6 7.644658 1.955342e+00
## 91 10.2 13.500812 -3.300812e+00
## 92 13.2 11.375487 1.824512e+00
## 93 13.6 14.389249 -7.892484e-01
## 94 12.2 10.438151 1.761849e+00
## 95 6.8 10.536242 -3.736242e+00
## 96 6.4 10.617671 -4.217670e+00
## 97 14.0 8.709231 5.290769e+00
## 98 6.0 6.001235 -1.235055e-03
## 99 11.2 10.644747 5.552533e-01
## 100 7.4 6.722280 6.777201e-01
## 101 8.2 13.654157 -5.454157e+00
## 102 13.4 14.017422 -6.174228e-01
## 103 17.6 13.053976 4.546025e+00
## 104 15.6 14.332258 1.267742e+00
## 105 4.2 6.042770 -1.842770e+00
## 106 14.0 10.682593 3.317407e+00
## 107 8.2 7.999778 2.002221e-01
## 108 8.2 12.552094 -4.352094e+00
## 109 21.8 14.292686 7.507313e+00
## 110 7.8 8.582041 -7.820404e-01
## 111 6.0 8.343737 -2.343737e+00
## 112 9.2 8.194452 1.005548e+00
## 113 6.4 7.582459 -1.182459e+00
## 114 17.4 14.403412 2.996588e+00
## 115 8.6 6.754897 1.845103e+00
## 116 5.8 7.167453 -1.367453e+00
## 117 8.0 7.723281 2.767194e-01
## 118 10.0 14.660207 -4.660207e+00
## 119 13.6 14.596500 -9.964993e-01
## 120 7.6 6.169289 1.430711e+00
## 121 8.2 7.880390 3.196094e-01
## 122 8.0 8.340124 -3.401239e-01
## 123 9.0 6.659345 2.340655e+00
## 124 23.4 14.215280 9.184720e+00
## 125 6.4 7.864002 -1.464002e+00
## 126 9.0 13.240818 -4.240818e+00
## 127 5.2 7.498141 -2.298141e+00
## 128 5.8 6.115133 -3.151327e-01
## 129 8.2 11.396106 -3.196107e+00
## 130 4.2 15.205634 -1.100563e+01
## 131 4.2 6.121779 -1.921779e+00
## 132 22.6 11.720766 1.087923e+01
## 133 6.8 6.093506 7.064944e-01
## 134 9.4 8.086294 1.313706e+00
## 135 8.0 14.693765 -6.693765e+00
## 136 11.0 11.349664 -3.496640e-01
## 137 11.8 11.534306 2.656946e-01
## 138 8.4 8.243185 1.568147e-01
## 139 15.8 12.998569 2.801431e+00
## 140 17.4 12.823034 4.576965e+00
## 141 6.0 8.445514 -2.445514e+00
## 142 4.4 6.028498 -1.628498e+00
## 143 12.6 14.578007 -1.978007e+00
## 144 10.6 14.636105 -4.036104e+00
## 145 9.2 7.534481 1.665519e+00
## 146 7.8 7.534481 2.655189e-01
## 147 6.0 8.273884 -2.273884e+00
## 148 7.2 6.131792 1.068207e+00
## 149 16.2 12.970539 3.229462e+00
## 150 5.2 6.424399 -1.224399e+00
## 151 6.6 7.598055 -9.980551e-01
## 152 6.0 6.276323 -2.763227e-01
## 153 6.2 8.580754 -2.380754e+00
## 154 13.6 11.488431 2.111569e+00
## 155 8.0 12.644147 -4.644147e+00
## 156 9.0 8.045075 9.549247e-01
## 157 11.2 12.893379 -1.693379e+00
## 158 10.2 14.565071 -4.365071e+00
## 159 6.2 6.007230 1.927699e-01
## 160 5.8 6.949630 -1.149630e+00
## 161 8.4 8.532963 -1.329634e-01
## 162 14.6 14.681413 -8.141237e-02
## 163 6.6 7.000992 -4.009924e-01
## 164 9.6 8.582447 1.017554e+00
## 165 5.2 7.797227 -2.597227e+00
## 166 10.4 8.593453 1.806547e+00
## 167 7.4 8.434236 -1.034236e+00
## 168 10.2 8.540487 1.659513e+00
## 169 6.8 6.011037 7.889633e-01
## 170 24.8 15.505213 9.294786e+00
## 171 19.0 14.595988 4.404012e+00
## 172 3.6 7.802669 -4.202670e+00
## 173 8.2 7.767130 4.328702e-01
## 174 7.0 6.856878 1.431218e-01
## 175 6.2 7.871803 -1.671803e+00
## 176 4.6 6.024248 -1.424248e+00
## 177 9.4 8.460099 9.399009e-01
## 178 6.6 8.643027 -2.043027e+00
## 179 5.6 7.477394 -1.877394e+00
## 180 6.4 6.077933 3.220671e-01
## 181 21.6 14.625048 6.974952e+00
## 182 15.2 15.505213 -3.052133e-01
## 183 12.4 14.267751 -1.867752e+00
## 184 5.8 6.038632 -2.386317e-01
## 185 13.6 14.681097 -1.081097e+00
## 186 19.4 14.567084 4.832915e+00
## 187 7.0 7.871803 -8.718032e-01
## 188 6.8 8.118681 -1.318681e+00
## 189 7.0 8.151121 -1.151121e+00
## 190 10.4 8.908466 1.491534e+00
## 191 7.0 7.782300 -7.823001e-01
## 192 14.2 14.711376 -5.113759e-01
## 193 8.8 8.053960 7.460404e-01
## 194 9.2 8.213257 9.867425e-01
## 195 6.4 8.005561 -1.605561e+00
## 196 6.8 6.132141 6.678593e-01
## 197 5.0 7.670591 -2.670591e+00
## 198 7.2 6.279358 9.206418e-01
## 199 7.6 9.626177 -2.026177e+00
## 200 8.2 8.314049 -1.140487e-01
## 201 8.6 6.978869 1.621131e+00
## 202 5.4 8.176950 -2.776949e+00
## 203 6.8 8.461626 -1.661626e+00
## 204 5.0 6.249640 -1.249640e+00
## 205 8.4 8.386907 1.309213e-02
## 206 24.8 14.636713 1.016329e+01
## 207 5.4 6.112184 -7.121835e-01
## 208 6.8 8.463866 -1.663865e+00
## 209 8.0 7.175164 8.248362e-01
## 210 6.2 8.505506 -2.305506e+00
## 211 6.4 8.152076 -1.752076e+00
## 212 7.0 8.295866 -1.295866e+00
## 213 6.6 9.598350 -2.998351e+00
## 214 7.2 8.485331 -1.285331e+00
## 215 9.0 9.748525 -7.485247e-01
## 216 7.0 9.336623 -2.336623e+00
## 217 7.8 9.559385 -1.759385e+00
## 218 11.2 14.699181 -3.499182e+00
## 219 12.2 8.100819 4.099181e+00
## 220 16.0 9.620612 6.379388e+00
## 221 11.0 14.353860 -3.353860e+00
## 222 13.4 14.524457 -1.124458e+00
## 223 6.0 9.448015 -3.448015e+00
## 224 8.2 14.391524 -6.191524e+00
## 225 18.2 14.434646 3.765355e+00
## 226 5.8 6.103081 -3.030804e-01
## 227 7.6 6.526578 1.073422e+00
## 228 6.6 7.347681 -7.476808e-01
## 229 7.0 6.829115 1.708847e-01
## 230 9.8 9.953834 -1.538336e-01
## 231 6.0 6.066767 -6.676699e-02
## 232 17.2 14.654966 2.545035e+00
## 233 9.8 8.582496 1.217504e+00
## 234 11.0 8.919558 2.080442e+00
## 235 6.2 8.914012 -2.714012e+00
## 236 12.0 14.672455 -2.672455e+00
## 237 10.4 8.551514 1.848486e+00
## 238 9.2 14.353860 -5.153860e+00
## 239 7.6 7.880390 -2.803905e-01
## 240 15.4 14.505711 8.942885e-01
## 241 6.0 6.141041 -1.410407e-01
## 242 8.4 8.227234 1.727653e-01
## 243 7.4 14.353860 -6.953860e+00
## 244 19.0 11.067380 7.932620e+00
## 245 7.4 9.487000 -2.087000e+00
## 246 9.6 14.686853 -5.086853e+00
## 247 16.2 13.623244 2.576756e+00
## 248 12.2 11.829915 3.700851e-01
## 249 10.8 13.933974 -3.133973e+00
## 250 4.8 6.075657 -1.275657e+00
## 251 14.0 10.601409 3.398591e+00
## 252 9.2 9.197428 2.571357e-03
## 253 5.4 6.003207 -6.032073e-01
## 254 5.6 7.060729 -1.460729e+00
## 255 7.4 9.275367 -1.875366e+00
## 256 5.0 8.210723 -3.210723e+00
## 257 10.0 9.258662 7.413377e-01
## 258 7.2 7.749697 -5.496974e-01
## 259 7.8 8.104119 -3.041192e-01
## 260 9.2 13.367232 -4.167232e+00
## 261 7.2 14.488184 -7.288185e+00
## 262 5.4 6.020570 -6.205699e-01
## 263 7.8 7.841898 -4.189804e-02
## 264 6.0 6.355638 -3.556380e-01
## 265 12.0 8.275928 3.724072e+00
## 266 8.6 7.107064 1.492936e+00
## 267 10.4 10.438151 -3.815096e-02
## 268 12.6 12.120473 4.795271e-01
## 269 8.0 11.062098 -3.062098e+00
## 270 8.6 7.464823 1.135177e+00
## 271 5.0 14.654152 -9.654152e+00
## 272 8.2 7.914493 2.855069e-01
## 273 9.2 7.323215 1.876785e+00
## 274 3.2 9.804079 -6.604079e+00
## 275 16.8 14.238196 2.561803e+00
## 276 6.0 6.002342 -2.342177e-03
## 277 12.2 14.514297 -2.314297e+00
## 278 9.0 14.435644 -5.435644e+00
## 279 16.2 14.595476 1.604525e+00
## 280 9.4 12.139451 -2.739451e+00
## 281 7.0 7.560589 -5.605886e-01
## 282 9.0 14.435145 -5.435145e+00
## 283 11.0 14.435644 -3.435644e+00
## 284 15.6 14.435644 1.164356e+00
## 285 9.0 9.581653 -5.816529e-01
## 286 15.4 11.072661 4.327339e+00
## 287 11.0 14.328492 -3.328492e+00
## 288 12.0 14.643621 -2.643621e+00
## 289 8.2 7.676862 5.231379e-01
## 290 6.4 8.379537 -1.979537e+00
## 291 29.2 14.632017 1.456798e+01
## 292 7.2 7.550160 -3.501600e-01
## 293 24.8 14.472600 1.032740e+01
## 294 10.6 9.631741 9.682589e-01
## 295 13.0 14.709716 -1.709716e+00
## 296 7.8 7.652395 1.476051e-01
## 297 13.0 14.666260 -1.666260e+00
## 298 10.0 14.528870 -4.528870e+00
## 299 7.0 6.986226 1.377433e-02
## 300 7.8 14.354062 -6.554062e+00
## 301 5.2 7.315368 -2.115369e+00
## 302 14.0 14.547008 -5.470076e-01
## 303 15.0 14.333695 6.663047e-01
## 304 4.0 8.259609 -4.259609e+00
## 305 11.8 7.511987 4.288013e+00
## 306 8.0 9.448015 -1.448015e+00
## 307 10.2 8.516921 1.683078e+00
## 308 14.2 8.275928 5.924072e+00
## 309 12.0 14.353517 -2.353517e+00
## 310 11.4 7.607623 3.792377e+00
## 311 12.0 12.848001 -8.480005e-01
## 312 6.6 8.207292 -1.607292e+00
## 313 8.2 7.323215 8.767846e-01
## 314 20.2 14.354348 5.845653e+00
## 315 14.6 14.402612 1.973883e-01
## 316 14.2 14.479729 -2.797292e-01
## 317 10.2 7.832254 2.367746e+00
## 318 10.0 14.596663 -4.596663e+00
## 319 20.2 14.354206 5.845794e+00
## 320 11.0 14.353860 -3.353860e+00
## 321 8.4 8.419673 -1.967304e-02
## 322 14.0 14.325163 -3.251633e-01
## 323 10.2 8.411002 1.788998e+00
## 324 9.8 8.508037 1.291963e+00
## 325 5.2 6.377820 -1.177820e+00
## 326 19.0 14.451422 4.548578e+00
## 327 6.2 8.248736 -2.048737e+00
## 328 8.2 8.360004 -1.600039e-01
## 329 6.6 14.199640 -7.599641e+00
## 330 7.2 7.363924 -1.639243e-01
## 331 6.4 9.626177 -3.226177e+00
## 332 15.2 11.665701 3.534299e+00
## 333 6.6 7.519837 -9.198367e-01
## 334 13.2 10.284826 2.915173e+00
## 335 27.0 14.957988 1.204201e+01
## 336 9.0 14.953536 -5.953536e+00
## 337 13.8 14.949104 -1.149104e+00
## 338 19.0 14.353778 4.646222e+00
## 339 9.4 14.498569 -5.098570e+00
## 340 9.8 8.450295 1.349705e+00
## 341 30.4 14.632932 1.576707e+01
## 342 8.4 7.598055 8.019447e-01
## 343 8.4 8.185584 2.144157e-01
## 344 14.2 8.654053 5.545947e+00
## 345 6.4 6.179211 2.207889e-01
## 346 9.2 7.971297 1.228703e+00
## 347 6.2 7.276204 -1.076205e+00
## 348 9.0 6.696727 2.303273e+00
## 349 7.0 14.475108 -7.475108e+00
## 350 8.2 6.085882 2.114118e+00
## 351 6.0 6.013272 -1.327184e-02
## 352 8.2 7.782208 4.177918e-01
## 353 8.2 8.324767 -1.247671e-01
## 354 9.8 14.353517 -4.553516e+00
## 355 10.0 14.353517 -4.353517e+00
## 356 6.2 6.377820 -1.778198e-01
## 357 7.6 8.692669 -1.092669e+00
## 358 5.8 9.859589 -4.059589e+00
## 359 11.4 14.325163 -2.925164e+00
## 360 6.8 7.633824 -8.338243e-01
## 361 17.0 14.333695 2.666305e+00
## 362 20.6 14.598005 6.001995e+00
## 363 12.4 12.485320 -8.532022e-02
## 364 7.6 6.304137 1.295863e+00
## 365 7.0 8.814277 -1.814277e+00
## 366 11.4 10.790328 6.096718e-01
## 367 9.2 8.159364 1.040635e+00
## 368 6.4 8.161945 -1.761945e+00
## 369 19.8 13.077452 6.722547e+00
## 370 8.8 8.527974 2.720258e-01
## 371 10.2 14.680461 -4.480462e+00
## 372 27.4 14.705254 1.269475e+01
## 373 5.0 14.127763 -9.127763e+00
## 374 21.0 14.402612 6.597388e+00
## 375 11.8 15.194758 -3.394758e+00
## 376 5.6 9.175168 -3.575168e+00
## 377 4.6 14.628789 -1.002879e+01
## 378 12.4 14.136265 -1.736265e+00
## 379 19.2 14.693614 4.506387e+00
## 380 8.0 14.320139 -6.320139e+00
## 381 7.6 8.557400 -9.573999e-01
## 382 8.8 8.478415 3.215851e-01
## 383 8.4 8.360004 3.999593e-02
## 384 8.2 6.760455 1.439544e+00
## 385 6.2 7.386105 -1.186105e+00
## 386 7.8 6.913525 8.864749e-01
## 387 7.6 14.488184 -6.888184e+00
## 388 20.2 14.353215 5.846786e+00
## 389 21.5 14.711502 6.788498e+00
## 390 18.2 14.258637 3.941363e+00
## 391 17.2 14.234422 2.965579e+00
## 392 6.0 6.410113 -4.101127e-01
## 393 8.4 8.084864 3.151360e-01
## 394 11.8 10.147269 1.652731e+00
## 395 6.0 8.548841 -2.548841e+00
## 396 4.2 9.113973 -4.913973e+00
## 397 7.4 8.869661 -1.469661e+00
## 398 6.2 7.000992 -8.009925e-01
## 399 7.6 8.583735 -9.837349e-01
## 400 11.2 13.338872 -2.138872e+00
## 401 7.0 6.141041 8.589593e-01
## 402 6.8 8.463866 -1.663865e+00
## 403 18.2 14.333695 3.866305e+00
## 404 6.4 8.584115 -2.184115e+00
## 405 8.0 7.566837 4.331632e-01
## 406 9.0 7.832254 1.167746e+00
## 407 8.8 14.415764 -5.615764e+00
## 408 9.2 11.503745 -2.303745e+00
## 409 8.0 8.580697 -5.806970e-01
## 410 5.8 6.020214 -2.202141e-01
## 411 20.4 13.926333 6.473667e+00
## 412 15.8 14.385521 1.414480e+00
## 413 18.0 9.637306 8.362694e+00
## 414 19.0 13.657218 5.342782e+00
## 415 9.2 13.702484 -4.502484e+00
## 416 16.0 13.772395 2.227605e+00
## 417 4.6 6.273301 -1.673301e+00
## 418 7.4 8.434236 -1.034236e+00
## 419 6.4 9.959372 -3.559372e+00
## 420 23.4 8.665082 1.473492e+01
## 421 6.2 6.301374 -1.013747e-01
## 422 17.8 14.711079 3.088921e+00
## 423 13.2 14.711079 -1.511079e+00
## 424 6.4 6.034712 3.652879e-01
## 425 6.2 6.001147 1.988532e-01
## 426 7.8 7.578447 2.215533e-01
## 427 8.8 8.532963 2.670372e-01
## 428 8.0 7.992833 7.166706e-03
## 429 8.0 7.856548 1.434522e-01
## 430 11.2 8.579158 2.620841e+00
## 431 6.6 6.147886 4.521137e-01
## 432 7.6 7.736804 -1.368040e-01
## 433 5.0 8.478073 -3.478073e+00
## 434 14.8 9.425737 5.374263e+00
## 435 4.0 7.451510 -3.451510e+00
## 436 6.0 6.001147 -1.146644e-03
## 437 6.2 8.314321 -2.114321e+00
## 438 13.8 14.685684 -8.856841e-01
## 439 36.8 14.595476 2.220452e+01
## 440 10.0 12.674452 -2.674452e+00
## 441 14.0 14.705098 -7.050978e-01
## 442 6.0 6.213853 -2.138529e-01
## 443 7.6 14.698099 -7.098099e+00
## 444 5.8 7.511987 -1.711987e+00
## 445 10.8 8.409508 2.390492e+00
## 446 26.2 14.705098 1.149490e+01
## 447 7.6 8.201292 -6.012918e-01
## 448 7.4 6.617179 7.828206e-01
## 449 5.8 6.114543 -3.145431e-01
## 450 9.2 14.705408 -5.505408e+00
## 451 5.6 6.002514 -4.025140e-01
## 452 10.4 9.637306 7.626939e-01
## 453 20.2 14.711079 5.488922e+00
## 454 11.4 11.976545 -5.765459e-01
## 455 9.0 14.354509 -5.354509e+00
## 456 7.0 6.129957 8.700431e-01
## 457 7.8 7.729189 7.081132e-02
## 458 10.4 6.320067 4.079933e+00
## 459 5.2 6.124243 -9.242430e-01
## 460 8.0 8.571676 -5.716759e-01
## 461 7.0 8.908466 -1.908466e+00
## 462 4.6 6.034080 -1.434080e+00
## 463 14.8 14.681097 1.189031e-01
## 464 9.8 8.541440 1.258560e+00
## 465 13.2 14.363987 -1.163987e+00
## 466 6.0 8.150418 -2.150418e+00
## 467 13.4 14.548435 -1.148435e+00
## 468 21.4 14.705254 6.694746e+00
## 469 8.4 13.460549 -5.060549e+00
## 470 16.4 14.471340 1.928660e+00
## 471 17.2 14.604603 2.595398e+00
## 472 7.2 8.415361 -1.215361e+00
## 473 12.0 10.075527 1.924473e+00
## 474 27.2 14.758165 1.244184e+01
## 475 5.4 7.560589 -2.160589e+00
## 476 5.8 8.254014 -2.454014e+00
## 477 6.2 6.021493 1.785069e-01
## 478 8.8 8.581965 2.180357e-01
## 479 7.0 6.274169 7.258310e-01
## 480 5.0 6.034712 -1.034712e+00
## 481 14.8 9.214126 5.585874e+00
## 482 6.0 8.107880 -2.107880e+00
## 483 12.8 10.312267 2.487733e+00
## 484 8.2 8.040917 1.590825e-01
## 485 16.0 14.420564 1.579436e+00
## 486 8.2 8.126693 7.330717e-02
## 487 8.8 8.061134 7.388660e-01
## 488 7.2 7.717273 -5.172730e-01
## 489 26.0 14.435644 1.156436e+01
## 490 7.0 7.909151 -9.091507e-01
## 491 7.4 7.221630 1.783703e-01
## 492 11.2 14.688281 -3.488281e+00
## 493 7.6 7.841898 -2.418983e-01
## 494 11.4 14.427251 -3.027251e+00
## 495 8.4 14.951318 -6.551319e+00
## 496 7.0 6.559619 4.403807e-01
## 497 6.2 13.409167 -7.209167e+00
## 498 7.8 6.087968 1.712033e+00
## 499 5.8 6.005137 -2.051373e-01
## 500 14.0 12.489797 1.510203e+00
## 501 5.8 6.117201 -3.172005e-01
## 502 9.6 8.235260 1.364740e+00
## 503 23.2 14.211397 8.988604e+00
## 504 8.2 6.323293 1.876707e+00
## 505 9.6 8.563296 1.036705e+00
## 506 6.6 14.353517 -7.753517e+00
## 507 10.6 12.275473 -1.675473e+00
## 508 9.0 8.559238 4.407617e-01
## 509 5.4 6.870904 -1.470904e+00
## 510 5.6 6.386898 -7.868983e-01
## 511 15.0 14.138379 8.616210e-01
## 512 10.2 8.045075 2.154924e+00
## 513 9.0 7.756636 1.243364e+00
## 514 12.7 14.578007 -1.878008e+00
## 515 13.2 10.470891 2.729109e+00
## 516 8.6 8.242551 3.574491e-01
## 517 6.4 6.010372 3.896278e-01
## 518 7.4 7.819936 -4.199362e-01
## 519 8.2 14.063853 -5.863853e+00
## 520 6.6 9.297640 -2.697640e+00
## 521 4.2 6.138791 -1.938791e+00
## 522 9.8 14.332165 -4.532165e+00
## 523 17.5 14.640118 2.859882e+00
## 524 8.4 11.104304 -2.704305e+00
## 525 8.2 8.265328 -6.532793e-02
## 526 7.0 8.323803 -1.323803e+00
## 527 8.2 7.930529 2.694705e-01
## 528 11.4 14.705098 -3.305098e+00
## 529 6.8 7.083384 -2.833834e-01
## 530 13.4 14.353860 -9.538605e-01
## 531 9.2 6.849901 2.350099e+00
## 532 14.4 14.323493 7.650699e-02
## 533 11.4 14.403250 -3.003251e+00
## 534 11.2 13.559777 -2.359777e+00
## 535 27.2 14.571624 1.262838e+01
## 536 8.4 8.460099 -6.009913e-02
## 537 9.0 8.324767 6.752330e-01
## 538 12.6 14.354348 -1.754347e+00
## 539 10.6 8.159364 2.440636e+00
## 540 8.6 14.389662 -5.789661e+00
## 541 7.2 6.677892 5.221074e-01
## 542 7.2 6.684139 5.158611e-01
## 543 5.3 14.711033 -9.411033e+00
## 544 6.6 7.602388 -1.002388e+00
## 545 16.0 7.936433 8.063567e+00
## 546 8.0 10.530803 -2.530803e+00
## 547 7.8 6.605417 1.194584e+00
## 548 8.8 13.866236 -5.066236e+00
## 549 18.2 14.402612 3.797389e+00
## 550 10.0 14.566961 -4.566961e+00
## 551 9.2 14.458528 -5.258528e+00
## 552 11.4 8.295866 3.104134e+00
## 553 8.0 8.825348 -8.253484e-01
## 554 16.0 14.389249 1.610751e+00
## 555 15.6 10.536242 5.063758e+00
## 556 7.0 6.870904 1.290964e-01
## 557 14.8 12.695981 2.104020e+00
## 558 4.0 13.711392 -9.711392e+00
## 559 24.0 14.354062 9.645938e+00
## 560 10.0 11.590104 -1.590104e+00
## 561 8.8 8.345158 4.548419e-01
## 562 6.8 7.770864 -9.708640e-01
## 563 7.0 7.574651 -5.746513e-01
## 564 5.4 8.183610 -2.783610e+00
## 565 5.4 7.030732 -1.630732e+00
## 566 23.6 12.029641 1.157036e+01
## 567 6.0 6.000636 -6.364998e-04
## 568 7.0 6.207772 7.922282e-01
## 569 36.6 13.530510 2.306949e+01
## 570 9.8 7.346778 2.453222e+00
## 571 27.2 14.333695 1.286631e+01
## 572 11.4 8.530594 2.869406e+00
## 573 6.8 7.782300 -9.822999e-01
## 574 10.0 8.074475 1.925525e+00
## 575 15.8 10.666381 5.133619e+00
## 576 9.0 10.822529 -1.822529e+00
## 577 5.4 6.024379 -6.243787e-01
## 578 33.4 14.389249 1.901075e+01
## 579 8.6 7.992920 6.070805e-01
## 580 9.0 7.804713 1.195287e+00
## 581 11.0 14.389249 -3.389249e+00
## 582 12.6 12.986582 -3.865817e-01
## 583 8.0 8.423936 -4.239362e-01
## 584 16.4 11.099035 5.300965e+00
## 585 10.2 7.868525 2.331474e+00
## 586 7.4 7.856548 -4.565477e-01
## 587 6.8 6.077363 7.226370e-01
## 588 8.0 7.938038 6.196163e-02
## 589 7.8 7.663157 1.368429e-01
## 590 9.2 14.333812 -5.133812e+00
## 591 6.8 7.142964 -3.429642e-01
## 592 24.0 14.602592 9.397408e+00
## 593 15.0 14.480047 5.199533e-01
## 594 6.0 8.445514 -2.445514e+00
## 595 7.8 8.281310 -4.813098e-01
## 596 7.0 10.295806 -3.295806e+00
## 597 8.4 7.914493 4.855067e-01
## 598 8.0 6.087933 1.912067e+00
## 599 6.0 7.982163 -1.982163e+00
## 600 6.4 14.374855 -7.974855e+00
## 601 7.6 6.309405 1.290595e+00
## 602 11.8 10.530803 1.269197e+00
## 603 5.6 6.244778 -6.447776e-01
## 604 5.2 6.182113 -9.821129e-01
## 605 7.0 7.496278 -4.962776e-01
## 606 7.7 8.027307 -3.273067e-01
## 607 8.8 8.003568 7.964324e-01
## 608 8.8 14.136265 -5.336265e+00
## 609 11.0 7.578447 3.421553e+00
## 610 11.4 14.574036 -3.174037e+00
## 611 6.0 6.048027 -4.802683e-02
## 612 5.6 8.094276 -2.494276e+00
## 613 8.4 7.298203 1.101797e+00
## 614 7.4 6.170179 1.229821e+00
## 615 11.2 10.125210 1.074790e+00
## 616 7.0 11.303041 -4.303041e+00
## 617 9.0 7.844457 1.155543e+00
## 618 7.0 7.832254 -8.322536e-01
## 619 8.6 7.819936 7.800641e-01
## 620 12.0 14.403250 -2.403250e+00
## 621 8.2 8.875203 -6.752033e-01
## 622 15.6 14.752867 8.471337e-01
## 623 6.2 6.842948 -6.429481e-01
## 624 5.2 6.202442 -1.002442e+00
## 625 18.6 10.328719 8.271281e+00
## 626 9.2 14.193701 -4.993701e+00
## 627 4.4 6.003984 -1.603984e+00
## 628 4.8 6.204043 -1.404043e+00
## 629 8.4 8.478073 -7.807360e-02
## 630 10.0 12.087141 -2.087141e+00
## 631 6.2 6.346970 -1.469707e-01
## 632 5.4 8.562728 -3.162728e+00
## 633 6.8 7.592862 -7.928623e-01
## 634 8.2 7.075818 1.124182e+00
## 635 9.2 7.807977 1.392023e+00
## 636 10.0 14.690668 -4.690668e+00
## 637 5.2 8.167359 -2.967360e+00
## 638 10.2 13.795019 -3.595020e+00
## 639 9.0 9.386749 -3.867492e-01
## 640 16.2 14.575293 1.624708e+00
## 641 9.8 8.580697 1.219303e+00
## 642 6.6 14.221069 -7.621070e+00
## 643 5.8 8.501835 -2.701835e+00
## 644 6.4 8.172775 -1.772775e+00
## 645 6.2 6.373324 -1.733246e-01
## 646 7.0 8.551514 -1.551514e+00
## 647 11.2 8.337189 2.862811e+00
## 648 7.2 6.251957 9.480433e-01
## 649 8.2 8.542812 -3.428126e-01
## 650 9.0 8.535658 4.643420e-01
## 651 4.4 6.836019 -2.436019e+00
## 652 4.6 6.044197 -1.444197e+00
## 653 8.2 8.522562 -3.225625e-01
## 654 7.2 6.026406 1.173593e+00
## 655 14.2 13.199856 1.000144e+00
## 656 8.4 8.323803 7.619681e-02
## 657 8.0 8.403997 -4.039974e-01
## 658 22.2 14.353778 7.846222e+00
## 659 7.6 8.386907 -7.869076e-01
## 660 5.4 7.386105 -1.986105e+00
## 661 6.4 8.435682 -2.035682e+00
## 662 17.6 13.419536 4.180465e+00
## 663 7.4 6.123504 1.276496e+00
## 664 6.8 6.400738 3.992624e-01
## 665 5.2 14.402612 -9.202612e+00
## 666 6.2 6.008500 1.915000e-01
## 667 5.2 6.229690 -1.029690e+00
## 668 9.8 8.337189 1.462812e+00
## 669 14.0 14.619493 -6.194932e-01
## 670 8.0 6.899230 1.100770e+00
## 671 16.8 9.937214 6.862785e+00
## 672 9.2 8.604464 5.955363e-01
## 673 8.0 10.114174 -2.114174e+00
## 674 17.6 11.473097 6.126903e+00
## 675 13.6 14.364237 -7.642367e-01
## 676 6.8 8.576348 -1.776348e+00
## 677 9.0 10.427227 -1.427227e+00
## 678 11.0 14.142593 -3.142593e+00
## 679 9.0 8.864120 1.358802e-01
## 680 14.0 13.244513 7.554873e-01
## 681 6.0 6.232496 -2.324958e-01
## 682 9.6 7.393976 2.206024e+00
## 683 9.2 12.458377 -3.258377e+00
## 684 7.2 9.653996 -2.453996e+00
## 685 11.6 8.309036 3.290964e+00
## 686 6.4 8.107880 -1.707880e+00
## 687 5.6 7.649569 -2.049570e+00
## 688 4.0 6.338420 -2.338420e+00
## 689 5.8 7.237187 -1.437187e+00
## 690 20.4 14.711350 5.688650e+00
## 691 19.8 14.595476 5.204523e+00
## 692 14.0 14.487640 -4.876399e-01
## 693 14.6 12.889275 1.710726e+00
## 694 10.2 8.319499 1.880500e+00
## 695 12.0 8.703709 3.296291e+00
## 696 7.0 7.247597 -2.475966e-01
## 697 6.2 6.321666 -1.216665e-01
## 698 13.6 14.676197 -1.076196e+00
## 699 18.0 13.081348 4.918652e+00
## 700 11.0 12.058447 -1.058447e+00
## 701 6.2 6.144903 5.509690e-02
## 702 5.0 6.002674 -1.002674e+00
## 703 10.4 7.386105 3.013894e+00
## 704 9.2 9.742967 -5.429671e-01
## 705 7.0 6.360016 6.399843e-01
## 706 5.6 7.411911 -1.811911e+00
## 707 15.2 15.495431 -2.954314e-01
## 708 14.4 14.330150 6.984981e-02
## 709 6.0 7.675851 -1.675851e+00
## 710 18.0 11.072661 6.927339e+00
## 711 17.0 14.705254 2.294746e+00
## 712 7.5 6.722280 7.777200e-01
## 713 12.8 14.353517 -1.553516e+00
## 714 7.0 8.498652 -1.498652e+00
## 715 5.4 8.053960 -2.653960e+00
## 716 4.6 6.143306 -1.543307e+00
## 717 9.2 7.827074 1.372926e+00
## 718 12.0 11.447491 5.525087e-01
## 719 19.2 13.976416 5.223584e+00
## 720 11.0 8.281370 2.718630e+00
## 721 10.2 6.159617 4.040382e+00
## 722 6.0 6.548482 -5.484817e-01
## 723 9.0 13.412629 -4.412629e+00
## 724 7.8 6.468877 1.331123e+00
## 725 13.8 14.681097 -8.810969e-01
## 726 8.2 14.333695 -6.133696e+00
## 727 8.6 6.424399 2.175602e+00
## 728 11.2 14.643621 -3.443621e+00
## 729 8.0 7.551191 4.488094e-01
## 730 13.0 7.550160 5.449840e+00
## 731 18.6 14.512872 4.087128e+00
## 732 8.2 6.690417 1.509583e+00
## 733 17.0 10.207865 6.792135e+00
## 734 8.2 6.323293 1.876707e+00
## 735 6.4 6.041366 3.586341e-01
## 736 12.4 14.403412 -2.003412e+00
## 737 6.6 8.579999 -1.979999e+00
## 738 15.0 14.389249 6.107512e-01
## 739 9.0 11.615367 -2.615367e+00
## 740 5.4 6.001378 -6.013779e-01
## 741 7.0 14.482624 -7.482624e+00
## 742 6.8 6.089796 7.102038e-01
## 743 14.5 14.490364 9.635703e-03
## 744 8.8 14.508464 -5.708463e+00
## 745 7.0 6.368859 6.311411e-01
## 746 6.2 8.521290 -2.321291e+00
## 747 10.6 7.508528 3.091472e+00
## 748 6.8 8.330544 -1.530544e+00
## 749 6.8 7.802669 -1.002669e+00
## 750 16.0 14.402612 1.597388e+00
## 751 7.2 6.877951 3.220487e-01
## 752 7.6 14.353778 -6.753779e+00
## 753 11.6 9.903960 1.696041e+00
## 754 11.0 14.640118 -3.640118e+00
## 755 8.0 7.331066 6.689338e-01
## 756 5.4 7.298203 -1.898203e+00
## 757 11.6 14.651396 -3.051396e+00
## 758 11.2 14.681727 -3.481727e+00
## 759 7.4 6.494646 9.053538e-01
## 760 6.4 6.131792 2.682076e-01
## 761 7.0 7.819637 -8.196367e-01
## 762 4.2 14.435145 -1.023514e+01
## 763 8.8 6.885022 1.914979e+00
## 764 6.2 8.577287 -2.377287e+00
## 765 6.6 6.605417 -5.416774e-03
## 766 8.4 8.485331 -8.533163e-02
## 767 4.0 13.270242 -9.270242e+00
## 768 7.0 6.099199 9.008008e-01
## 769 8.0 12.818858 -4.818858e+00
## 770 20.0 14.228726 5.771274e+00
## 771 19.0 14.354062 4.645938e+00
## 772 5.4 11.995895 -6.595895e+00
## 773 24.6 14.747607 9.852394e+00
## 774 5.0 6.025380 -1.025380e+00
## 775 11.0 9.865137 1.134863e+00
## 776 5.0 7.534532 -2.534532e+00
## 777 8.6 7.252774 1.347227e+00
## 778 12.3 14.752867 -2.452867e+00
## 779 6.0 6.009686 -9.685971e-03
## 780 10.6 10.025784 5.742166e-01
## 781 7.0 8.557579 -1.557579e+00
## 782 13.6 14.435644 -8.356435e-01
## 783 7.6 14.290933 -6.690933e+00
## 784 13.6 14.452459 -8.524582e-01
## 785 7.8 8.061134 -2.611340e-01
## 786 6.8 6.010372 7.896279e-01
## 787 5.4 7.504134 -2.104134e+00
## 788 10.0 10.213368 -2.133683e-01
## 789 5.6 6.364423 -7.644226e-01
## 790 5.6 8.698189 -3.098189e+00
## 791 28.0 13.737824 1.426218e+01
## 792 10.0 6.377820 3.622180e+00
## 793 9.4 8.265328 1.134672e+00
## 794 5.2 6.141041 -9.410409e-01
## 795 14.0 11.590104 2.409896e+00
## 796 9.8 14.487640 -4.687640e+00
## 797 6.4 8.266359 -1.866359e+00
## 798 7.4 6.255588 1.144412e+00
## 799 6.0 6.293459 -2.934588e-01
## 800 5.0 15.197470 -1.019747e+01
## 801 8.4 6.316854 2.083145e+00
## 802 17.4 14.097409 3.302590e+00
## 803 9.4 9.080610 3.193900e-01
## 804 9.2 14.203577 -5.003577e+00
## 805 17.8 13.553312 4.246688e+00
## 806 7.0 8.880746 -1.880746e+00
## 807 11.0 13.850113 -2.850113e+00
## 808 6.4 7.949235 -1.549235e+00
## 809 11.0 11.365164 -3.651644e-01
## 810 7.0 6.096391 9.036089e-01
## 811 10.2 10.003658 1.963421e-01
## 812 9.2 8.891832 3.081679e-01
## 813 6.2 6.781508 -5.815079e-01
## 814 8.0 6.956911 1.043089e+00
## 815 5.8 6.453780 -6.537794e-01
## 816 24.0 10.972012 1.302799e+01
## 817 8.0 7.943562 5.643772e-02
## 818 7.8 8.185584 -3.855838e-01
## 819 20.8 14.711033 6.088966e+00
## 820 8.6 12.708849 -4.108849e+00
## 821 11.6 11.780466 -1.804659e-01
## 822 6.6 6.599582 4.175813e-04
## 823 8.6 6.499892 2.100109e+00
## 824 12.0 14.680461 -2.680461e+00
## 825 5.8 6.559619 -7.596191e-01
## 826 7.8 7.744405 5.559473e-02
## 827 10.4 10.345161 5.483832e-02
## 828 14.2 14.353860 -1.538603e-01
## 829 5.0 8.797673 -3.797673e+00
## 830 7.8 13.014489 -5.214489e+00
## 831 7.0 8.423936 -1.423936e+00
## 832 7.0 7.770864 -7.708642e-01
## 833 8.0 6.405411 1.594589e+00
## 834 6.8 8.350189 -1.550189e+00
## 835 19.6 14.504016 5.095984e+00
## 836 18.0 14.389249 3.610751e+00
## 837 14.2 14.217214 -1.721447e-02
## 838 15.2 14.458528 7.414719e-01
## 839 5.0 7.578447 -2.578447e+00
## 840 17.6 14.354062 3.245938e+00
## 841 15.2 14.708154 4.918459e-01
## 842 24.8 15.495431 9.304568e+00
## 843 8.8 14.501863 -5.701862e+00
## 844 6.0 7.941229 -1.941229e+00
## 845 8.0 9.659559 -1.659559e+00
## 846 6.0 6.892115 -8.921148e-01
## 847 9.0 6.329985 2.670015e+00
## 848 5.6 6.010372 -4.103724e-01
## 849 9.4 9.158475 2.415245e-01
## 850 7.0 14.470076 -7.470076e+00
## 851 8.0 8.581434 -5.814341e-01
## 852 7.4 8.483558 -1.083558e+00
## 853 16.4 12.556516 3.843484e+00
## 854 12.0 8.352229 3.647771e+00
## 855 5.8 6.029823 -2.298225e-01
## 856 14.8 9.141784 5.658216e+00
## 857 4.0 6.136558 -2.136558e+00
## 858 8.0 6.592513 1.407487e+00
## 859 4.4 6.065796 -1.665796e+00
## 860 29.6 13.613864 1.598614e+01
## 861 11.2 9.492569 1.707431e+00
## 862 11.0 8.324767 2.675233e+00
## 863 4.2 13.273899 -9.073899e+00
## 864 6.8 8.753429 -1.953429e+00
## 865 13.8 14.562237 -7.622365e-01
## 866 14.0 14.389249 -3.892488e-01
## 867 5.2 6.005137 -8.051377e-01
## 868 5.4 8.273884 -2.873884e+00
## 869 7.4 7.717273 -3.172727e-01
## 870 14.4 14.632017 -2.320174e-01
## 871 8.2 8.024538 1.754621e-01
## 872 6.2 6.006673 1.933270e-01
## 873 7.0 7.691643 -6.916433e-01
## 874 22.6 14.588557 8.011443e+00
## 875 12.8 14.354348 -1.554348e+00
## 876 7.7 6.020570 1.679430e+00
## 877 6.6 6.510474 8.952560e-02
## 878 7.0 7.782208 -7.822080e-01
## 879 5.2 6.285470 -1.085470e+00
## 880 5.8 6.611282 -8.112820e-01
## 881 7.0 6.044197 9.558027e-01
## 882 5.9 6.000121 -1.001207e-01
## 883 9.8 8.295866 1.504134e+00
## 884 7.0 10.741921 -3.741921e+00
## 885 8.4 8.309036 9.096316e-02
## 886 4.4 7.776160 -3.376160e+00
## 887 22.0 14.353778 7.646222e+00
## 888 5.2 14.222990 -9.022991e+00
## 889 6.8 6.000060 7.999406e-01
## 890 4.2 6.001147 -1.801147e+00
## 891 9.0 8.714753 2.852471e-01
## 892 8.2 14.613558 -6.413558e+00ggplot(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)
Natural 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-16ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black")+
stat_smooth(method = "lm",
formula = y~ns(x, knots = knot_value_manual_3),
lty = 1,se=F)
Smoothing 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()
Pada Smoothing Splin dilakukan smoothing di hampir semua bagian.
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
Dapat dilihat bahwa, jika lambda semakin kecil maka plot semakin kompleks sehingga menghasilkan hubungan x dan y yang pemulusannya lebih banyak, semakin besar lambda maka ujung plot akan terlihat seperti garis terlihat seperti linear.
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.20355Didapatkan nilai GCV sebesar 14.20355, semakin kecil nilai GCV maka semakin bagus.
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()
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
Lakukan analisis regresi non linier pada data berikut:
library(ISLR)
AutoData = Auto %>% select(mpg, horsepower, origin)
tibble(AutoData)## # A tibble: 392 x 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 rowsX <- AutoData$horsepower
Y <- AutoData$mpgRegresi Linier
lin.mod <-lm( Y~X)
plot(X,Y)
lines(X,lin.mod$fitted.values,col="red")
Polynomial
pol.mod <- lm( Y~X+I(X^2))
ix <- sort(X,index.return=T)$ix
plot(X,Y)
lines(X[ix], pol.mod$fitted.values[ix],col="blue", cex=2)
summary(pol.mod)##
## Call:
## lm(formula = Y ~ X + I(X^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.7135 -2.5943 -0.0859 2.2868 15.8961
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.9000997 1.8004268 31.60 <2e-16 ***
## X -0.4661896 0.0311246 -14.98 <2e-16 ***
## I(X^2) 0.0012305 0.0001221 10.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.374 on 389 degrees of freedom
## Multiple R-squared: 0.6876, Adjusted R-squared: 0.686
## F-statistic: 428 on 2 and 389 DF, p-value: < 2.2e-16Fungsi Tangga
#regresi fungsi tangga
range(X)## [1] 46 230step.mod = lm(Y ~ cut(X,5))
summary(step.mod)##
## Call:
## lm(formula = Y ~ cut(X, 5))
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.3033 -3.0220 -0.5413 2.4074 16.6394
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.3033 0.4272 73.27 <2e-16 ***
## cut(X, 5)(82.8,120] -8.2813 0.5642 -14.68 <2e-16 ***
## cut(X, 5)(120,156] -15.2427 0.7210 -21.14 <2e-16 ***
## cut(X, 5)(156,193] -17.5922 1.0036 -17.53 <2e-16 ***
## cut(X, 5)(193,230] -18.5340 1.3767 -13.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.719 on 387 degrees of freedom
## Multiple R-squared: 0.6382, Adjusted R-squared: 0.6345
## F-statistic: 170.7 on 4 and 387 DF, p-value: < 2.2e-16ggplot(AutoData,aes(x= X, y=Y)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = Y~cut(X,5),
lty = 1, col = "blue",se = F)+
theme_bw()## Warning: 'newdata' had 80 rows but variables found have 392 rows## Warning: Computation failed in `stat_smooth()`:
## arguments imply differing number of rows: 80, 392
Komparasi 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 2363.324
## 2 Poly (ordo=2) 2274.354
## 3 Tangga (breaks=3) 2335.820Model yang terbaik dapat dilihat dari nilai AIC yang terkecil. Berdasarkan output komperasi model di atas, dapat dilihat bahwa model terbaik adalah model Polynomial dengan ordo=2 dengan nilai AIC sebesar 2274.354.