Analisis Regresi Berganda

Secara sederhana dikatakan analisis regresi berganda karena pada umumnya memiliki lebih dari satu peubah bebas atau peubah independent/prediktor.

Load Package

Karena data yang digunakan terismpan dalam format xlsx/ file excel.

#memanggil package
library(readxl)

Input data

#input data
dat=read_excel("E:\\Praktikum-5.xlsx",sheet="Sheet1")
dat
## # A tibble: 25 x 4
##       NO     y    x1    x2
##    <dbl> <dbl> <dbl> <dbl>
##  1     1  16.7     7   560
##  2     2  11.5     3   220
##  3     3  12.0     3   340
##  4     4  14.9     4    80
##  5     5  13.8     6   150
##  6     6  18.1     7   330
##  7     7   8       2   110
##  8     8  17.8     7   210
##  9     9  79.2    30  1460
## 10    10  21.5     5   605
## # ... with 15 more rows

Plot data

#melihat plot
plot(dat[,-1]) 

#x1 vs y
plot(dat$x1,dat$y,main="Jumlah Minuman (X1) dan Lama Pengantaran (Y)",
     xlab="Jumlah Minuman (X1)",ylab="Lama Pengantaran (Y)")

#x2 vs y
plot(dat$x2,dat$y,main="Jarak Tempuh (X2) dan Lama Pengantaran (Y)",
     xlab="Jarak Tempuh (X2)",ylab="Lama Pengantaran (Y)")

Pendugaan dengan MKT

#membuat model regresi
reg_ganda=lm(y~x1+x2,data=dat)
coef(reg_ganda)
## (Intercept)          x1          x2 
##  2.34123115  1.61590721  0.01438483
n=nrow(dat)
p=length(coef(reg_ganda))
k=p-1
X=cbind(rep(1,n),dat$x1,dat$x2)
y=dat$y

# solve(A) --> invers(A)
# t(A) --> transpose(A)
beta_duga = solve(t(X) %*% X) %*% t(X) %*% y;beta_duga
##            [,1]
## [1,] 2.34123115
## [2,] 1.61590721
## [3,] 0.01438483

Membuat tabel sidik ragam (ANOVA)

#derajat bebas
db_r<-k
db_s<-n-(k+1)
y_duga<-reg_ganda$fitted.values
y_rataan<-mean(dat$y)
# Manual ANOVA
JKR=sum((y_duga-y_rataan)^2)
JKR
## [1] 5550.811
JKS=sum((y-y_duga)^2)
JKS
## [1] 233.7317
JKT=sum((y-y_rataan)^2)
JKT
## [1] 5784.543
KTR=JKR/db_r
KTR
## [1] 2775.405
KTS=JKS/db_s
KTS
## [1] 10.62417
tabel_anova=data.frame(SK=c("Regresi","Sisaan","Total"),
                       db=c(db_r,db_s,db_r+db_s),
                       JK=c(JKR,JKS,JKT),
                       KT=c(KTR,KTS,""),
                       Fhit=c(KTR/KTS,"",""))
tabel_anova
##        SK db        JK               KT             Fhit
## 1 Regresi  2 5550.8109 2775.40546128972 261.235108660564
## 2  Sisaan 22  233.7317 10.6241671554797                 
## 3   Total 24 5784.5426
#mengeluarkan anova melalui fungsi anova
anova(reg_ganda)
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## x1         1 5382.4  5382.4 506.619 < 2.2e-16 ***
## x2         1  168.4   168.4  15.851 0.0006312 ***
## Residuals 22  233.7    10.6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uji Simultan

#Uji simultan (Uji F)
Fhit=KTR/KTS
Fhit
## [1] 261.2351
Ftabel=qf(0.95,db_r,db_s)
Ftabel
## [1] 3.443357

Uji Parsial (tiap-tiap beta)

#simpangan baku sisaan (akar KTS,atau dugaan bagi sigma)
akar_KTS=sqrt(KTS)
akar_KTS
## [1] 3.259473
#standar error koefisien (se)
diag_XX=diag(solve(t(X)%*%X))
se_beta=as.vector(akar_KTS) * sqrt(diag_XX)
se_beta
## [1] 1.096730168 0.170734918 0.003613086
#uji parsial (tiap-tiap beta)
thitung=beta_duga/se_beta
thitung
##          [,1]
## [1,] 2.134738
## [2,] 9.464421
## [3,] 3.981313
ttabel=qt(0.975,db_s)
ttabel
## [1] 2.073873
#uji parsial B0
thitung[1,]
## [1] 2.134738
if (thitung[1,] >  ttabel){
  print("tolak H0")
  }else {
  print("Terima H0")}
## [1] "tolak H0"
#uji parsial B1
thitung[2,]
## [1] 9.464421
if (thitung[2,] >  ttabel){
  print("tolak H0")
}else {
  print("Terima H0")}
## [1] "tolak H0"
#uji parsial B2
thitung[3,]
## [1] 3.981313
if (thitung[3,] >  ttabel){
  print("tolak H0")
}else {
  print("Terima H0")}
## [1] "tolak H0"

Uji Kebaikan Model

#ukuran kebaikan model
R_sq=JKR/JKT
R_sq
## [1] 0.9595937
R_sq_adj=1-(KTS/(JKT/(n-1)))
R_sq_adj
## [1] 0.9559205

Selang Kepercayaan bagi parameter

#SK bagi beta
confint(reg_ganda,level=0.95)
##                   2.5 %     97.5 %
## (Intercept) 0.066751987 4.61571030
## x1          1.261824662 1.96998976
## x2          0.006891745 0.02187791
#manual SK
BB_beta=beta_duga-qt(0.975,22)*se_beta
BA_beta=beta_duga+qt(0.975,22)*se_beta
tabel_SK_beta=data.frame(Parameter=c("Beta0","Beta1","Beta2"),
                         BB=BB_beta,BA=BA_beta)
tabel_SK_beta
##   Parameter          BB         BA
## 1     Beta0 0.066751987 4.61571030
## 2     Beta1 1.261824662 1.96998976
## 3     Beta2 0.006891745 0.02187791

Selang Kepercayaan Individu dan Rataan

#amatan baru
baru=data.frame(x1=11, x2=90)

#SK untuk individu y
predict.lm(reg_ganda,baru,interval ="prediction" )
##        fit      lwr      upr
## 1 21.41084 13.86086 28.96083
#Sk untuk rataan y
predict.lm(reg_ganda,baru,interval ="confidence" )
##        fit      lwr      upr
## 1 21.41084 18.04806 24.77363
#manual Sk individu dan rataan y

x0=matrix(c(1,11,90),3,1)
x0
##      [,1]
## [1,]    1
## [2,]   11
## [3,]   90
y0_duga=t(x0) %*%  beta_duga
y0_duga
##          [,1]
## [1,] 21.41084
a=t(x0) %*% solve(t(X)%*%X) %*% x0
a
##           [,1]
## [1,] 0.2474787
BB_individu_y=y0_duga - ttabel*(akar_KTS*sqrt(1+a))
BA_individu_y=y0_duga + ttabel*(akar_KTS*sqrt(1+a))

#manual Sk individu dan rataan y
BB_rataan_y=y0_duga - ttabel*(akar_KTS*sqrt(a))
BA_rataan_y=y0_duga + ttabel*(akar_KTS*sqrt(a))


sk_y=data.frame(sk=c("individu","rataan"),
                lower=c(BB_individu_y,BB_rataan_y),
                upper=c(BA_individu_y,BA_rataan_y)) 
sk_y
##         sk    lower    upper
## 1 individu 13.86086 28.96083
## 2   rataan 18.04806 24.77363