Model Reduksi Hasil Korelasi

model2 = lm(formula = Y ~ X1+X2+X4+X6+X7+X8+X9+X12, data = dt)
summary(model2)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X4 + X6 + X7 + X8 + X9 + X12, data = dt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4006  -1.4476   0.5167   2.3660   6.7547 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  5.806e+01  5.405e+01   1.074   0.3023  
## X1          -7.717e-01  6.726e-01  -1.147   0.2719  
## X2           9.612e-01  3.126e+00   0.308   0.7633  
## X4           4.282e-01  6.137e-01   0.698   0.4976  
## X6          -1.058e-01  1.554e-01  -0.681   0.5078  
## X7          -2.572e+00  1.415e+00  -1.817   0.0923 .
## X8           9.913e-04  5.710e-04   1.736   0.1062  
## X9          -5.652e-05  7.422e-05  -0.762   0.4599  
## X12          1.318e-01  2.054e-01   0.642   0.5321  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.906 on 13 degrees of freedom
## Multiple R-squared:  0.6744, Adjusted R-squared:  0.4741 
## F-statistic: 3.366 on 8 and 13 DF,  p-value: 0.02559

Diperoleh model sebagai berikut
\[ \hat Y = 0.5806 - 0.7717X_1 - 0.9612X_2 + 0.4282X_4 - 0.1058X_6 - 2.572X_7 + 0.0009X_8 - 0.00005X_9 + 0.1318X_{12} + e \]

Model Hasil Reduksi Multikolinearitas

model3 = lm(formula = Y ~ X1+X2+X4+X6+X7, data = dt)
summary(model3)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X4 + X6 + X7, data = dt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5808  -2.2473  -0.8028   2.5361   8.3073 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  63.6722    25.7074   2.477   0.0248 *
## X1           -0.8045     0.5253  -1.531   0.1452  
## X2            2.9462     2.6568   1.109   0.2838  
## X4            0.7075     0.5500   1.286   0.2167  
## X6           -0.1301     0.1351  -0.963   0.3498  
## X7           -2.2446     1.3801  -1.626   0.1234  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.974 on 16 degrees of freedom
## Multiple R-squared:  0.588,  Adjusted R-squared:  0.4593 
## F-statistic: 4.568 on 5 and 16 DF,  p-value: 0.008863

Diperoleh model sebagai berikut
\[ \hat Y = 63.6722 - 0.8045X_1 + 2.9462X_2 + 0.7075X_4 - 0.1301X_6 - 2.2446X_7 + e \]

1. Nilai Harapan Sisaan sama dengan Nol

H0: Nilai harapan sisaan sama dengan nol
H1: Nilai harapan sisaan tidak sama dengan nol

t.test(model3$residuals,mu = 0,conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  model3$residuals
## t = -7.6223e-17, df = 21, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.924963  1.924963
## sample estimates:
##     mean of x 
## -7.055446e-17

Uji t.tes tersebut menunjukkan hasil p-value = 1 > alpha = 0.05, maka tak tolak H0, nilai harapan sisaan sama dengan nol pada taraf nyata 5%. Asumsi terpenuhi.

2. Ragam Sisaan Homogen

\(H0:var[e]=sigma2I\) (ragam sisan homogen)
\(H1:var[e] != sigma2I\) (ragam siaan tidak homogen)

bptest(model3)

## 
##  studentized Breusch-Pagan test
## 
## data:  model3
## BP = 3.7978, df = 5, p-value = 0.5789

Uji ini sering disebut dengan uji homokesdasitas yang dilakukan dengan yang dilakukan dengan uji Breusch-Pagan. Karena p-value = 0.5789 > alpha = 0.05, maka tak tolak H0, ragam sisaan homogen pada taraf nyata 5%. Asumsi terpenuhi.

3. Sisaan Saling Bebas

\(H0:E[ei,ej]=0\) (sisaan saling bebas/tidak ada autokorelasi)
\(H1:E[ei,ej] != 0\) (sisaan tidak saling bebas/ada autokorelasi)

dwtest(model3)

## 
##  Durbin-Watson test
## 
## data:  model3
## DW = 2.2648, p-value = 0.634
## alternative hypothesis: true autocorrelation is greater than 0

Uji ini sering disebut dengan uji autokorelasi yang dilakukan dengan Durbin watson. Karena p-value = 0.634 (pada DW test) > alpha = 0.05, maka tak tolak H0, sisaan saling bebas pada taraf nyata 5%, sehingga asumsi terpenuhi.

Internally Studentized Residuals Plot (ri)

ggplot(hasil) +
  geom_point(aes(y = `ri`, x = `index`),
             color="blue",size=2) +
  ylab("ri") +
  xlab("index") +
  ggtitle("Plot ri") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))

Cook Distance Plot (ci)

ggplot(hasil) +
  geom_point(aes(y = `ci`, x = `index`),
             color="blue",size=2) +
  ylab("ci") +
  xlab("index") +
  ggtitle("Plot ci") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))

Difference in Fits Plot (DFFITSi)

ggplot(hasil) +
  geom_point(aes(y = `DFFITSi`, x = `index`),
             color="blue",size=2) +
  ylab("DFFITSi") +
  xlab("index") +
  ggtitle("Plot DFFITSi") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))

Hadi’s Influence Measure (Hi)

ggplot(hasil) +
  geom_point(aes(y = `hadi`, x = `index`),
             color="blue",size=2) +
  ylab("Hi") +
  xlab("index") +
  ggtitle("Plot Hi") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))

Potential Residual Plot (P-R Plot)

ggplot(hasil) +
  geom_point(aes(y = `potential`, x = `residual`),
             color="blue",size=2) +
  ylab("potential") +
  xlab("residual") +
  ggtitle("P-R Plot") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))

Pencilan atau Outlier

for (i in 1:dim(hasil)[1]){
  absri <- abs(hasil$ri)
  pencilan <- which(absri > 2)
}
pencilan

## [1] 21

Pencilan adalah titik amatan yang menjauh pada tren Y walau X nya masih berada pada selang. Dalam data tersebut terdapat dua titik observasi yang terdeteksi sebagai sebuah pencilan, yaitu amatan ke-21 dengan ri bernilai -2.7629 yang absolutnya lebih besar dari 2..

Titik Leverage

titik_leverage <- vector("list", dim(hasil)[1])
for (i in 1:dim(hasil)[1]) {
  cutoff <- 2 * 6 / 22
  titik_leverage[[i]] <- which(hii > cutoff)
}
leverages <- unlist(titik_leverage)
titik_leverage <- sort(unique(leverages))
titik_leverage

## [1] 22

Titik leverage adalah nilai amatan yang X nya jauh lebih besar dari X lainnya dan terletak hampir dalam garis regresi tetapi juga menjauh. Terdapat dua titik observasi yang termasuk sebagai titik leverage pada data tersebut, yaitu amatan ke-22 dengan nilai hii 0.9051. Kedua titik ini memiliki nilai hii yang lebih besar dari 2*p/n = 0.5454, dengan p = 6 adalah banyaknya parameter dan n = 22 adalah banyaknya amatan.

Amatan Berpengaruh

amatan_berpengaruh <- vector("list", dim(hasil)[1])
for (i in 1:dim(hasil)[1]) {
  cutoff <- 2 * sqrt((6 / 22))
  amatan_berpengaruh[[i]] <- which(abs(DFFITSi) > cutoff)
}
berpengaruh <- unlist(amatan_berpengaruh)
amatan_berpengaruh <- sort(unique(berpengaruh))
amatan_berpengaruh

## [1] 10 21 22

Amatan berpengaruh adalah nilai amatan yang memengaruhi hasil dugaan parameter regresi, R-square, dan uji hipotesis jika disisihkan. Hasil pada data tersebut menunjukkan bahwa amatan ke-10, 21, dan 22 termasuk sebagai amatan berpengaruh yang memiliki nilai DFFITSi -1.0587, -1.2123, dan -1.1085. Absolut dari nilai tersebut lebih besar dari 2*sqrt(p/n) = 1.0445.

#gambaran amatan berpengaruh dengan pencilan dan laverage
ols_plot_resid_lev(model3)

data <- dt[-c(21, 22),]

model_berpengaruh <- lm(Y~X1+X2+X4+X6+X7, data=data)

summary(model3)$r.squared

## [1] 0.5880444

summary(model_berpengaruh)$r.squared

## [1] 0.6563153

Berdasarkan hasil tersebut, nilai adj-Rsquare lebih kecil ketika pencilan dan leverage dikeluarkan dari data dibandingkan yang tidak. Oleh karenanya kedua titik tersebut disisihkan.

#pembentukan data baru tanpa amatan 9 dan 19
dt2 <- subset(dt, select = -c(X3, X5, X8, X9, X10, X11, X12))
dt2 <- dt2[-c(21, 22),]
#pemodelan rlb terbaru
model4 <- lm(Y~.,data= dt2)
summary(model4)

## 
## Call:
## lm(formula = Y ~ ., data = dt2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6156 -2.8722  0.1116  2.3695  7.4795 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 61.30733   39.37661   1.557   0.1418  
## X1          -1.03474    0.50058  -2.067   0.0577 .
## X2           4.72896    3.61451   1.308   0.2118  
## X4           1.14999    0.60974   1.886   0.0802 .
## X6          -0.09898    0.12744  -0.777   0.4503  
## X7          -2.34888    1.29452  -1.814   0.0911 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.322 on 14 degrees of freedom
## Multiple R-squared:  0.6563, Adjusted R-squared:  0.5336 
## F-statistic: 5.347 on 5 and 14 DF,  p-value: 0.005904

anova(model4)

## Analysis of Variance Table
## 
## Response: Y
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## X1         1 185.445 185.445  9.9263 0.007083 **
## X2         1  55.202  55.202  2.9548 0.107647   
## X4         1  87.575  87.575  4.6876 0.048143 * 
## X6         1 109.740 109.740  5.8740 0.029501 * 
## X7         1  61.508  61.508  3.2923 0.091085 . 
## Residuals 14 261.552  18.682                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

bs <- ols_step_best_subset(model4)
bs

##    Best Subsets Regression   
## -----------------------------
## Model Index    Predictors
## -----------------------------
##      1         X6             
##      2         X4 X7          
##      3         X1 X4 X7       
##      4         X1 X2 X4 X7    
##      5         X1 X2 X4 X6 X7 
## -----------------------------
## 
##                                                    Subsets Regression Summary                                                    
## ---------------------------------------------------------------------------------------------------------------------------------
##                        Adj.        Pred                                                                                           
## Model    R-Square    R-Square    R-Square     C(p)       AIC        SBIC        SBC         MSEP        FPE       HSP       APC  
## ---------------------------------------------------------------------------------------------------------------------------------
##   1        0.4100      0.3772      0.2939    8.0337    124.9836    67.4983    127.9708    499.1876    27.4392    1.4673    0.7211 
##   2        0.5203      0.4639      0.3734    5.5394    122.8431    66.2259    126.8260    431.2049    24.6940    1.3421    0.6490 
##   3        0.6097      0.5365      0.4462    3.8983    120.7187    65.8920    125.6973    374.2412    22.2763    1.2376    0.5854 
##   4        0.6415      0.5459      0.3664    4.6033    121.0193    67.6842    126.9937    368.3101    22.7352    1.2992    0.5975 
##   5        0.6563      0.5336      0.2921    6.0000    122.1755    70.1935    129.1457    380.2561    24.2870    1.4371    0.6383 
## ---------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria 
##  SBIC: Sawa's Bayesian Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
##  MSEP: Estimated error of prediction, assuming multivariate normality 
##  FPE: Final Prediction Error 
##  HSP: Hocking's Sp 
##  APC: Amemiya Prediction Criteria

null<-lm(Y~ 1, data=dt2) # 1 here means the intercept 
full<-lm(Y~., data=dt2)
 
step(full, scope=list(lower=null, upper=full),data=datanew, direction='backward', trace=0)

## 
## Call:
## lm(formula = Y ~ X1 + X4 + X7, data = dt2)
## 
## Coefficients:
## (Intercept)           X1           X4           X7  
##     87.2720      -0.9279       0.6748      -2.8960

step(null, scope=list(lower=null, upper=full),data=datanew, direction='forward', trace=0)

## 
## Call:
## lm(formula = Y ~ X6 + X1 + X7 + X4 + X2, data = dt2)
## 
## Coefficients:
## (Intercept)           X6           X1           X7           X4           X2  
##    61.30733     -0.09898     -1.03474     -2.34888      1.14999      4.72896

step(null, scope=list(lower=null, upper=full),data=datanew, direction='both', trace=0)

## 
## Call:
## lm(formula = Y ~ X1 + X7 + X4, data = dt2)
## 
## Coefficients:
## (Intercept)           X1           X7           X4  
##     87.2720      -0.9279      -2.8960       0.6748

modell1 <-  lm(Y~X1+X4+X7, data = dt2)
modell2 <- lm(Y~X1+X2+X4+X6+X7, data = dt2)
summary(modell1)$adj.r.squared

## [1] 0.5365348

summary(modell2)$adj.r.squared

## [1] 0.5335708

ols_aic(modell1)

## [1] 120.7187

ols_aic(modell2)

## [1] 122.1755

model <- modell1
model

## 
## Call:
## lm(formula = Y ~ X1 + X4 + X7, data = dt2)
## 
## Coefficients:
## (Intercept)           X1           X4           X7  
##     87.2720      -0.9279       0.6748      -2.8960