Nama kelompok:

Faradilla ilena (8309141002)

Qatrunnada Azkia (18309141007)
Hanan Albarr (18309141024)
Reynaldi Komtua Naibaho
Prodi : Statistika G 2018

A. Data mathdata
> mathdata <- read.table("D:\\data\\mathproficiencydata.txt",header = T)
> str(mathdata)
'data.frame': 40 obs. of 7 variables:
$ State : Factor w/ 40 levels "Alabama","Arizona",..: 1 2 3 4 5 6 8 7 9 10 ...
$ Mathprof: int 252 259 256 256 267 270 261 231 255 258 ...
$ Parents : int 75 75 77 78 78 79 75 47 75 73 ...
$ Homelib : int 78 73 77 68 85 86 83 76 73 80 ...
$ Reading : int 34 41 28 42 38 43 32 24 31 36 ...
$ TVwatch : int 18 12 20 11 9 12 18 33 19 17 ...
$ Absences: int 18 26 23 28 25 22 28 37 27 22 ...

Apakah ada oulier?

> library(MASS)
> cooks <- cooks.distance(ols)
> plot(cooks,ylab = "Cooks distances")

Terlihat ada outlier disekitar data 35 dan 40

> data.frame(State=mathdata$State,Cook.distance=round(cooks,6))
> library(MASS)
> cook.dist<-cooks.distance(ols)
> r<-stdres(ols)
> a<-cbind(mathdata,cook.dist,r)
> n<-nrow(mathdata)
> a[cook.dist>-4/n,]
State Mathprof Parents Homelib Reading TVwatch Absences cook.dist r
1 Alabama 252 75 78 34 18 18 9.487414e-05 -0.05463857
2 Arizona 259 75 73 41 12 26 6.088905e-03 0.40076248
3 Arkansas 256 77 77 28 20 23 6.255947e-02 1.39338192
4 California 256 78 68 42 11 28 7.298034e-04 0.10336573
5 Colorado 267 78 85 38 9 25 5.877966e-03 -0.57826419
6 Connecticut 270 79 86 43 12 22 2.313419e-05 0.03170589
7 Delaware 261 75 83 32 18 28 9.312647e-03 0.64984635
8 D.C. 231 47 76 24 33 37 7.191353e-01 1.39145359
9 Florida 255 75 73 31 19 27 3.487150e-02 1.44485143
10 Georgia 258 73 80 36 17 22 3.247754e-03 0.48432464
 11 Guam 231 81 64 32 20 28 5.688087e-01 -2.57511006
12 Hawaii 251 78 69 36 23 26 1.732292e-01 1.63516146
13 Idaho 272 84 84 48 7 21 2.157953e-02 -0.81378050
14 Illinois 260 78 82 43 14 21 2.196530e-02 -0.84876406
15 Indiana 267 81 84 37 11 23 1.777606e-04 -0.14909914
16 Iowa 278 83 88 43 8 20 2.131173e-03 0.36961066
17 Kentucky 256 79 78 36 14 23 2.069387e-03 -0.56122535
18 Louisiana 246 73 76 36 19 27 1.227221e-02 -1.05099520
19 Maryland 260 75 83 34 19 27 5.713268e-03 0.51136186
20 Michigan 264 77 84 31 14 25 4.357220e-03 0.43041654
21 Minnesota 276 83 88 36 7 20 2.846258e-03 0.31747066
22 Montana 280 83 88 44 6 21 1.765449e-03 0.32178103
23 Nebraska 276 85 88 42 9 19 1.806907e-04 0.09634409
24 New_Hampshire 273 83 88 40 7 22 5.877759e-03 -0.64905351
25 New_Jersey 269 79 84 41 13 23 2.603513e-03 0.44592687
26 New_Mexico 256 77 72 40 11 27 5.004062e-03 -0.35669699
27 New_York 261 76 79 35 17 29 1.034525e-02 0.82447217
28 North_Carolina 250 74 78 37 21 25 3.838767e-03 -0.36283139
29 North_Dakota 281 85 90 41 6 14 1.115378e-02 0.53466186
30 Ohio 264 79 84 36 11 22 2.274054e-03 -0.47521241
31 Oklahoma 263 78 78 37 14 22 5.365197e-03 0.83894760
32 Oregon 271 81 82 41 9 31 1.643958e-04 0.06038777
33 Pennsylvania 266 80 86 34 10 24 9.485219e-03 -0.62175058
34 Rhode_Island 260 78 80 38 12 28 7.052856e-03 -0.71406068
35 Texas 258 77 70 34 15 18 3.291897e-01 2.12282260
36 Virgin_Islands 218 63 76 23 27 22 1.205460e+00 -3.91733625
37 Virginia 264 78 82 33 16 24 9.440962e-03 0.90701530
38 West_Virginia 256 82 80 36 16 25 2.005284e-02 -0.91229850
39 Wisconsin 274 81 86 38 8 21 2.308753e-03 0.39381906
40 Wyoming 272 85 86 43 7 23 1.191973e-02 -0.91307933
Outlier pengamatan ke j dikatakan berpengaruh apabila D i > 4/n dengan n adalah banyaknya
pengamatan. Dari hasil R di atas dapat diketahui bahwa terdapat outlier yang berpengaruh adalah :
D.C, Guam, Hawaii, Texas, Virgin Island.

Karena ada outlier maka ada remidial. Sehingga dilakukan dengan robust regression. Dengan manual
R

Fungsi bobot

> ols <- lm(mathdata$Mathprof~mathdata$Parents+mathdata$Homelib+mathdata$Reading+mathdata$TVwatch+mathdata$Absences,data =

mathdata)
> ai <- (abs(ols$residuals-median(ols$residuals)))
> ai
1 2 3 4 5 6 7 8 9 10
0.73322417 1.43615532 6.25007136 0.01083135 3.36678993 0.31287137 2.74774223 3.61011214 6.78711754 1.98219576
11 12 13 14 15 16 17 18 19 20
11.49207350 6.84025886 4.39026890 4.58049421 1.23670743 1.39249009 3.36938449 5.83031331 2.06351507 1.65317864
21 22 23 24 25 26 27 28 29 30
1.07710600 1.14513591 0.01083135 3.75395318 1.79254875 2.15965369 3.68822268 2.23280247 2.06598665 2.90050680
31 32 33 34 35 36 37 38 39 40
3.85243131 0.18720312 3.52746383 4.08412224 8.85531322 17.48257509 4.15227693 4.96170526 1.51833800 5.08564460
> ei <- ols$residuals
> MAD <- median(abs(ols$residuals-median(ols$residuals)))/0.6745
> MAD
[1] 4.645883
> ui <- ei/MAD
> abs(ui)
1 2 3 4 5 6 7 8 9 10 11
0.05677825 0.41016835 1.44633642 0.09871270 0.62363831 0.03370030 0.69247998 0.87810020 1.56193256 0.52770045 2.37255942
12 13 14 15 16 17 18 19 20 21 22
1.57337092 0.84393635 0.88488126 0.16515016 0.40076966 0.62419678 1.15389772 0.54520397 0.45688139 0.33288503 0.34752808
23 24 25 26 27 28 29 30 31 32 33
0.10337547 0.70697300 0.48688001 0.36380910 0.89491305 0.37955396 0.54573597 0.52327352 0.93025803 0.06074968 0.65822246
34 35 36 37 38 39 40
0.77804002 2.00709992 3.66198084 0.99479809 0.96693477 0.42785771 0.99361201
> head(ui <- ols$residuals/MAD)
> head(abs(ui))
1 2 3 4 5 6
0.05677825 0.41016835 1.44633642 0.09871270 0.62363831 0.03370030

 Dengan huber

Tabel

No Ei Ui wi
1 -0.2637851 -0.05677825 1
2 1.9055943 0.41016835 1
3 6.7195104 1.44633642 0.929935
-
-
40 -4.6162056 -0.99361201 1
−1
Lalu hasil yang di cari di atas di masukan ke rumus b =( W WX ) X t WY untuk menentukan
¿ t

estimator yang baru

> rr.huber <-
rlm(mathdata$Mathprof~mathdata$Parents+mathdata$Homelib+mathdata$Reading+mathdata$TVwatch+m
athdata$Absences,data = mathdata,psi=psi.huber)
> summary(rr.huber)

Call: rlm(formula = mathdata$Mathprof ~ mathdata$Parents + mathdata$Homelib +

mathdata$Reading + mathdata$TVwatch + mathdata$Absences,
data = mathdata, psi = psi.huber)
Residuals:
Min 1Q Median 3Q Max
-22.1969 -3.1202 0.9977 2.5250 5.8669

Coefficients:
Value Std. Error t value
(Intercept) 178.1624 28.3043 6.2945
mathdata$Parents 0.3573 0.2008 1.7795
mathdata$Homelib 0.7703 0.1404 5.4875
mathdata$Reading 0.1495 0.2101 0.7118
mathdata$TVwatch -0.8218 0.2754 -2.9846
mathdata$Absences -0.0121 0.2059 -0.0588

Residual standard error: 4.094 on 34 degrees of freedom

> plot(rr.huber$w,ylab="Huber Weight")

 Dengan bisquare
Tabel

No ei Ui Wi
1 -0.2637851 -0.05677825 0.999
2 1.9055943 0.41016835 0.984
3 6.7195104 1.44633642 0.818
-
-
40 -4.6162056 -0.99361201 0.91
−1
Lalu hasil yang di cari di atas di masukan ke rumus b =( W WX ) X t WY untuk menentukan
¿ t

estimator yang baru

> rr.bsq <-
rlm(mathdata$Mathprof~mathdata$Parents+mathdata$Homelib+mathdata$Reading+mathdata$TVwatch+m
athdata$Absences,data = mathdata,psi = psi.bisquare)
> summary(rr.bsq)

Call: rlm(formula = mathdata$Mathprof ~ mathdata$Parents + mathdata$Homelib +

mathdata$Reading + mathdata$TVwatch + mathdata$Absences,
data = mathdata, psi = psi.bisquare)
Residuals:
Min 1Q Median 3Q Max
-24.9476 -3.2096 0.5676 2.2121 4.0740

Coefficients:
Value Std. Error t value
(Intercept) 186.5746 25.0714 7.4417
mathdata$Parents 0.4344 0.1779 2.4424
mathdata$Homelib 0.6734 0.1243 5.4162
mathdata$Reading 0.0116 0.1861 0.0621
 mathdata$TVwatch -0.7559 0.2439 -3.0991
mathdata$Absences -0.0924 0.1824 -0.5066

Residual standard error: 3.518 on 34 degrees of freedom

> plot(rr.bsq$w,ylab="Huber Weight")

Kesimpulan
Pada data mathdata menghasilakan residual standard eror dari ols =5.628, huber=4.094, bisquare =
3.518. dari hasil tersebut dapat disimpulkan bahwa regresi robust dengan bisquare paling bagus
digunakan karena residuals standart erornya paling kecil.

B. Data duncan
> library(car)
> str(Duncan)
'data.frame': 45 obs. of 4 variables:
$ type : Factor w/ 3 levels "bc","prof","wc": 2 2 2 2 2 2 2 2 3 2 ...
$ income : int 62 72 75 55 64 21 64 80 67 72 ...
$ education: int 86 76 92 90 86 84 93 100 87 86 ...
$ prestige : int 82 83 90 76 90 87 93 90 52 88 ...
> modell <- lm(Duncan$prestige ~ Duncan$income+Duncan$education,data=Duncan)
> summary(modell)

Call:
lm(formula = Duncan$prestige ~ Duncan$income + Duncan$education,
data = Duncan)

Residuals:
Min 1Q Median 3Q Max
-29.538 -6.417 0.655 6.605 34.641

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.06466 4.27194 -1.420 0.163
Duncan$income 0.59873 0.11967 5.003 1.05e-05 ***
Duncan$education 0.54583 0.09825 5.555 1.73e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.37 on 42 degrees of freedom

Multiple R-squared: 0.8282, Adjusted R-squared: 0.82
 F-statistic: 101.2 on 2 and 42 DF, p-value: < 2.2e-

> library(MASS)
> cooks <- cooks.distance(modell)
> plot(cooks,ylab = "Cooks distances")

Terlihat adanya oulier

Fungsi bobot
> ai <- (abs(modell$residuals-median(modell$residuals)))
> ai
accountant pilot architect author chemist
3.3469100 3.8179209 0.2883799 0.6452958 10.1494444
minister professor dentist reporter engineer
33.9866235 9.3286070 7.0719555 30.1925880 3.3595818
undertaker lawyer physician welfare.worker teacher
3.1284266 4.5853564 3.9604775 5.9880329 0.0000000
conductor contractor factory.owner store.manager banker
20.6519862 25.1146957 19.9193930 1.2465907 5.9505205
bookkeeper mail.carrier insurance.agent store.clerk carpenter
12.2532321 19.3499793 25.2744516 23.2448861 13.2824920
electrician RR.engineer machinist auto.repairman plumber
8.9820961 8.6293532 23.3889945 6.2295931 5.5800307
gas.stn.attendant coal.miner streetcar.motorman taxi.driver truck.driver
9.4001145 12.3980941 14.9283990 0.3493785 2.3508367
machine.operator barber bartender shoe.shiner cook
5.9199937 1.6386544 12.4530134 6.2577107 1.0194557
soda.clerk watchman janitor policeman waiter
12.1497500 7.4142445 1.6977467 0.3989515 6.8464865
> ei <- modell$residuals
> MAD <- median(abs(modell$residuals-median(modell$residuals)))/0.6745
> MAD
[1] 10.15046
> ui <- ei/MAD
> ui
> abs(ui)
accountant pilot architect author chemist
0.3942196776 0.4406225791 0.0929003686 0.0009167949 1.0643896158
minister professor dentist reporter engineer
3.4127733262 0.9835226186 0.6322228354 2.9100140289 0.3954680725
undertaker lawyer physician welfare.worker teacher
0.2437155005 0.3872488477 0.4546669270 0.5254372910 0.0644898475
conductor contractor factory.owner store.manager banker
1.9700960240 2.5387315341 2.0269023330 0.1873010747 0.6507213565
bookkeeper mail.carrier insurance.agent store.clerk carpenter
1.1426702161 1.8418253176 2.4254906145 2.2255425238 1.3730502147
electrician RR.engineer machinist auto.repairman plumber
0.9493851643 0.9146337475 2.3687194455 0.6782149401 0.4852418597
gas.stn.attendant coal.miner streetcar.motorman taxi.driver truck.driver
0.8615876741 1.2859213761 1.4062214614 0.0300698880 0.1671091429
machine.operator barber bartender shoe.shiner cook
0.6477139278 0.2259262858 1.1623522080 0.5520053217 0.1649242613
soda.clerk watchman janitor policeman waiter
1.1324753932 0.6659443592 0.1027682323 0.1037936272 0.6100101525
> head(abs(ui))
accountant pilot architect author chemist minister
0.3942196776 0.4406225791 0.0929003686 0.0009167949 1.0643896158 3.4127733262

 Dengan huber
 No ei ui wi
1 4.001511779 0.3942196776 0.98589
2 4.472522658 0.440622579 0,982387
3 0.942981643 0.092900368 0,999213
4 0.009305892 0.000916794 0.99999992
5 10.804046136 1.064389615 0.89943283
6 34.641225280 3.412773326 0.22030373
.
.
.
39 1.053553248 0.103793627 0.9990186
40 -6.191884751 -0.610010152 0.9663807
−1
Lalu hasil yang di cari di atas di masukan ke rumus b ¿=( W t WX ) X WY untuk menentukan
t

estimator yang baru

Dengan R
> hb <- rlm(Duncan$prestige ~Duncan$income+Duncan$education,data=Duncan, psi=
psi.huber)
> summary(hb)

Call: rlm(formula = Duncan$prestige ~ Duncan$income + Duncan$education,

data = Duncan, psi = psi.huber)
Residuals:
Min 1Q Median 3Q Max
-30.120 -6.889 1.291 4.592 38.603

Coefficients:
Value Std. Error t value
(Intercept) -7.1107 3.8813 -1.8320
Duncan$income 0.7014 0.1087 6.4516
Duncan$education 0.4854 0.0893 5.4380

Residual standard error: 9.892 on 42 degrees of freedom

> plot(hb$w,ylab="Huber Weight")

 Dengan bisquare

No ei ui wi
1 4.001511779 0.3942196776 0.98589
2 4.472522658 0.440622579 0,982387
3 0.942981643 0.092900368 0,999213
 4 0.009305892 0.000916794 0.99999992
5 10.804046136 1.064389615 0.89943283
6 34.641225280 3.412773326 0.22030373
.
.
.
39 1.053553248 0.103793627 0.9990186
40 -6.191884751 -0.610010152 0.9663807
−1
Lalu hasil yang di cari di atas di masukan ke rumus b ¿=( W t WX ) X WY untuk menentukan
t

estimator yang baru

Dengan R
> bsq <- rlm(Duncan$prestige ~Duncan$income+Duncan$education,data=Duncan, psi=
psi.bisquare)
> summary(bsq)

Call: rlm(formula = Duncan$prestige ~ Duncan$income + Duncan$education,

data = Duncan, psi = psi.bisquare)
Residuals:
Min 1Q Median 3Q Max
-29.949 -6.580 1.509 4.421 42.657

Coefficients:
Value Std. Error t value
(Intercept) -7.4121 3.8770 -1.9118
Duncan$income 0.7902 0.1086 7.2761
Duncan$education 0.4186 0.0892 4.6942

Residual standard error: 9.558 on 42 degrees of freedom

>plot(bsq$w,ylab="Huber Weight")

Kesimpulan
Pada data mathdata menghasilakan residual standard eror dari modell(ols) =13.37, huber=9.892,
bisquare = 9.558. dari hasil tersebut dapat disimpulkan bahwa regresi robust dengan bisquare paling
bagus digunakan karena residuals standart erornya paling kecil.

C. Data mesin

Banyak item yang rusak dihasilkan oleh mesin(Y ) diketahui linear dengan pengaturan kecepatan
mesin ( X ). Data berikut dikumpulkn dari catatan control kualitas terbaru.

i 1 2 3 4 5 6 7 8 9 10 11 12
 x 200 400 300 400 200 300 300 400 200 400 200 300
y 28 75 37 53 22 58 40 96 46 52 30 69
Lakukan perhitungan manual regeresi robust

Dengan huber
> x <- c(200,400,300,400,200,300,300,400,200,400,200,300)
> y <- c(28,75,37,53,22,58,40,96,46,52,30,69)
> data_baru<- data.frame(y,x)
> mode <- lm(y~x,data_baru)
> #interaksi 1
> ytopi <- -5.7500+0.1875*x
> ei <- y-ytopi
> MAD <- (1/0.6745)*median(abs(e-median(e)))
> ui <- e/MAD
> ui
[1] -0.2408929 0.3693690 -0.8672143 -1.0438690 -0.6263214 0.4817857 -0.6745000 1.7183690
0.9153929 -1.1081071
[11] -0.1124167 1.1884048
> wh <- c(1,1,1,1,1,1,1,1.345/1.718369,1,1,1,1)
> wh1<-cbind(wh)
> xm<-cbind(1,x)
> b <- solve(t(xm)%*%solve(t(diag(c(wh1))))%*%xm)%*%(t(xm)%*%(diag(c(wh1))))%*%y
> b
[,1]
5.3016189
x 0.1401359
>
> #interaksi 2
> ytopi2 <- 5.3016189+0.1401359*x
> e2 <- y-ytopi2
> MAD <- (1/0.6745)*median(abs(e2-median(e2)))
> u2 <- e2/MAD
> u2
[1] -0.5523879 1.4143510 -1.0721010 -0.8661880 -1.1743531 1.1047772 -0.7611184 3.5912292
1.3135077 -0.9698489
[11] -0.3450662 2.2450467
> wh2 <- c(1,1.345/(abs(1.4143510)),1,1,1,1,1,1.345/abs((3.5912292)),1,1,1,1.345/(abs(2.2450467)))
> wh22<-cbind(wh2)
> xm<-cbind(1,x)
> b2 <- solve(t(xm)%*%solve(t(diag(c(wh22))))%*%xm)%*%(t(xm)%*%(diag(c(wh22))))%*%y
> b2
[,1]
27.65843979
x 0.02596945
>
> #interaksi 3
> ytopi3 <- 27.65843979 +0.02596945*x
> e3 <- y-ytopi3
> MAD <- (1/0.6745)*median(abs(e3-median(e3)))
> u3 <- e3/MAD
> u3
[1] -0.23046087 1.75511576 0.07365152 0.71022816 -0.51543022 1.07104423 0.21613620
2.75250847 0.62444716
[10] 0.66273327 -0.13547109 1.59348804
> wh3<-c(1,1.345/(abs(1.677308760)),1,1,1,1,1,1.345/(abs(2.75250847)),1,1,1,1.345/
(abs(1.59348804)))
> wh33<-cbind(wh3)
> xm<-cbind(1,x)
> b3 <- solve(t(xm)%*%solve(t(diag(c(wh22))))%*%xm)%*%(t(xm)%*%(diag(c(wh22))))%*%y
> b3
 [,1]
27.65843979
x 0.02596945

> cbind(ei,ui,wh1,e2,u2,wh22,e3,u3,wh33)
ei ui wh e2 u2 wh2 e3 u3 wh3
[1,] -3.75 -0.2408929 1.000000 -5.328799 -0.5523879 1.0000000 -4.852330 -0.23046087 1.0000000
[2,] 5.75 0.3693690 1.000000 13.644021 1.4143510 0.9509662 36.953780 1.75511576 0.8018798
[3,] -13.50 -0.8672143 1.000000 -10.342389 -1.0721010 1.0000000 1.550725 0.07365152 1.0000000
[4,] -16.25 -1.0438690 1.000000 -8.355979 -0.8661880 1.0000000 14.953780 0.71022816 1.0000000
[5,] -9.75 -0.6263214 1.000000 -11.328799 -1.1743531 1.0000000 -10.852330 -0.51543022 1.0000000
[6,] 7.50 0.4817857 1.000000 10.657611 1.1047772 1.0000000 22.550725 1.07104423 1.0000000
[7,] -10.50 -0.6745000 1.000000 -7.342389 -0.7611184 1.0000000 4.550725 0.21613620 1.0000000
[8,] 26.75 1.7183690 0.782719 34.644021 3.5912292 0.3745236 57.953780 2.75250847 0.4886452
[9,] 14.25 0.9153929 1.000000 12.671201 1.3135077 1.0000000 13.147670 0.62444716 1.0000000
[10,] -17.25 -1.1081071 1.000000 -9.355979 -0.9698489 1.0000000 13.953780 0.66273327 1.0000000
[11,] -1.75 -0.1124167 1.000000 -3.328799 -0.3450662 1.0000000 -2.852330 -0.13547109 1.0000000
[12,] 18.50 1.1884048 1.000000 21.657611 2.2450467 0.5990967 33.550725 1.59348804 0.8440603

> hb <- rlm(y~x,dat,psi = psi.huber)

> plot(hb$w,ylab="Huber Weight")

Dengan bisquare
> x <- c(200,400,300,400,200,300,300,400,200,400,200,300)
> y <- c(28,75,37,53,22,58,40,96,46,52,30,69)
> data_baru<- data.frame(y,x)
> mode <- lm(y~x,data_baru)
> #interaksi 1
> ytopi <- -5.7500+0.1875*x
> ei <- y-ytopi
> MAD <- (1/0.6745)*median(abs(e-median(e)))
> ui <- e/MAD
> ui
[1] -0.2408929 0.3693690 -0.8672143 -1.0438690 -0.6263214 0.4817857 -0.6745000 1.7183690
0.9153929 -1.1081071
[11] -0.1124167 1.1884048
> wh <- lapply(ui, function(x)
+ if(abs(x)<=4.685) {(1-(x/4.685)^2)^2}
+ else {0})
> wh1<-cbind(wh)
> xm<-cbind(1,x)
> b <- solve(t(xm)%*%solve(t(diag(c(wh1))))%*%xm)%*%(t(xm)%*%(diag(c(wh1))))%*%y
> b
[,1]
7.7983591
x 0.1143934
> #interaksi 2
> ytopi2 <- 7.7983591+0.1143934*x
> e2 <- y-ytopi2
> MAD <- (1/0.6745)*median(abs(e2-median(e2)))
> u2 <- e2/MAD
> u2
[1] -0.28751083 2.30309036 -0.54949305 -0.05968357 -0.93190372 1.70588206 -0.22729661
4.55846547 1.64566784
[10] -0.16708239 -0.07271320 2.88726903
> wh2 <- lapply(u2, function(x)
+ if(abs(x)<=4.685) {(1-(x/4.685)^2)^2}
+ else {0})
> wh22<-cbind(wh2)
> xm<-cbind(1,x)
> b2 <- solve(t(xm)%*%solve(t(diag(c(wh22))))%*%xm)%*%(t(xm)%*%(diag(c(wh22))))%*%y
> b2
[,1]
62.6280523
x -0.1553515
> #interaksi 3
> ytopi3 <- 62.6280523 -0.1553515*x
> e3 <- y-ytopi3
> MAD <- (1/0.6745)*median(abs(e3-median(e3)))
> u3 <- e3/MAD
> u3
[1] -0.1213960 2.5424825 0.7157810 1.7918087 -0.3261252 1.4323332 0.8181456 3.2590347
0.4927917 1.7576872
[11] -0.0531529 1.8076701
> wh3 <- lapply(u3, function(x)
+ if(abs(x)<=4.685) {(1-(x/4.685)^2)^2}
+ else {0})
> wh33<-cbind(wh3)
> xm<-cbind(1,x)
> b3 <- solve(t(xm)%*%solve(t(diag(c(wh22))))%*%xm)%*%(t(xm)%*%(diag(c(wh22))))%*%y
> b3
[,1]
62.6280523
x -0.1553515
> cbind(ei,ui,wh1,e2,u2,wh22,e3,u3,wh33)
ei ui wh e2 u2 wh2 e3 u3 wh3
[1,] -3.75 -0.2408929 0.9947194 -2.677039 -0.2875108 0.992482 -3.557752 -0.121396 0.9986576
[2,] 5.75 0.369369 0.9876069 21.44428 2.30309 0.5750812 74.51255 2.542482 0.4977193
[3,] -13.5 -0.8672143 0.9326467 -5.116379 -0.5494931 0.9726764 20.9774 0.715781 0.9538605
[4,] -16.25 -1.043869 0.9031752 -0.5557191 -0.05968357 0.9996754 52.51255 1.791809 0.7288498
[5,] -9.75 -0.6263214 0.9645752 -8.677039 -0.9319037 0.9224333 -9.557752 -0.3261252 0.9903322
[6,] 7.5 0.4817857 0.9789614 15.88362 1.705882 0.752417 41.9774 1.432333 0.8217979
[7,] -10.5 -0.6745 0.9589748 -2.116379 -0.2272966 0.995298 23.9774 0.8181456 0.9399381
[8,] 26.75 1.718369 0.7490412 42.44428 4.558465 0.002839549 95.51255 3.259035 0.2663556
[9,] 14.25 0.9153929 0.9251045 15.32296 1.645668 0.7684525 14.44225 0.4927917 0.9779946
[10,] -17.25 -1.108107 0.891244 -1.555719 -0.1670824 0.9974579 51.51255 1.757687 0.7383019
[11,] -1.75 -0.1124167 0.9988488 -0.6770391 -0.0727132 0.9995183 -1.557752 -0.0531529 0.9997426
[12,] 18.5 1.188405 0.8754517 26.88362 2.887269 0.3846477 52.9774 1.80767 0.7244153

> bsq <- rlm(y~x,dat,psi= psi.bisquare)

> plot(bsq$w,ylab = "huber weight")