Tugas 2 - Analisis Regresi

1. Input Data

library(readr)
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, union

# membaca file csv
dtraw <- read.csv("C:/Users/LENOVO/Downloads/data liver.csv", header = TRUE )
head(dtraw)

##   No    X1    X2   X3    X4    X5   X6      Y
## 1  1 16.36  8.90 3.47  6.02 57.42 1.11 158.76
## 2  2 26.68 21.22 3.53 12.07 61.38 1.36 197.19
## 3  3 12.49 16.62 2.00  8.88 67.42 1.47 144.73
## 4  4  8.45 22.86 6.71  7.46 69.94 1.31 140.06
## 5  5 10.19 14.23 4.75  2.06 65.68 1.25 129.71
## 6  6 19.53 17.35 1.95  7.54 59.63 1.14 162.59

2. Membuat Data Frame

# membuat subset data df dan mengubah nama kolomnya
df <-  dtraw[, c("Y", "X1", "X2", "X4"), drop = FALSE]
colnames(df) <- c("y", "x1", "x2", "x4")
df <- as.data.frame(df)
print(df)

##         y    x1    x2    x4
## 1  158.76 16.36  8.90  6.02
## 2  197.19 26.68 21.22 12.07
## 3  144.73 12.49 16.62  8.88
## 4  140.06  8.45 22.86  7.46
## 5  129.71 10.19 14.23  2.06
## 6  162.59 19.53 17.35  7.54
## 7  178.48 20.65 10.48  4.88
## 8  120.90 22.96 14.23  3.69
## 9  191.24 21.22 21.64 11.94
## 10 150.03  8.11  3.16  8.82
## 11 173.44 24.74  7.84  3.68
## 12 211.98 11.38 15.71  7.20
## 13 193.49 15.82 15.04  9.89
## 14 164.04  8.36  9.01  3.40
## 15 156.97 12.04  9.72  6.03
## 16 208.36 10.97  4.58  5.55
## 17 154.62  7.97  9.33  4.17
## 18 137.38  7.46  6.11  2.99
## 19 180.15 29.09 15.71  9.35
## 20 228.47 10.30  8.54 10.78
## 21 153.62  7.82  4.41  4.19
## 22 121.31 14.71  6.29  6.16
## 23 157.37  8.54  6.73  5.52
## 24 211.27 23.05 11.34  3.00
## 25 178.16 13.12  5.86 10.92
## 26 174.89  7.41  9.11  5.50
## 27 142.98 14.59  5.59  3.75
## 28 165.59  8.52  6.52  6.92
## 29 141.54 18.97  6.35  5.61
## 30 238.22 35.41 36.36 15.00
## 31 138.42  4.55  1.27  2.83
## 32 247.45 22.59 28.70 10.35
## 33 140.27  9.21  4.55  7.92
## 34 216.06 18.32 11.61  8.07
## 35 144.18  5.69  6.88  2.78
## 36 156.22 11.21 11.92 10.29

# cek tipe data df (perlu karena ada error sebelumnya)
class(df)

## [1] "data.frame"

# identifikasi banyak n amatan dan p parameter untuk perhitungan manual
n <- 36 #banyak amatan
p <- 4 #banyak parameter(k+1)

# menyematkan df untuk mempermudah pemanggilan kolom 
attach(df)

3. Eksplorasi Data

# melihat plot antar peubah
par(mfrow=c(1,3))
plot(x1, y, main='Scatterplot X1 vs Y', pch=19, col="#F65100")
plot(x2, y, main='Scatterplot X2 vs Y', pch=19, col="#0144aa")
plot(x4, y, main='Scatterplot X4 vs Y', pch=19, col="#262626")

Ketiga plot eksplorasi memberikan gambaran kecenderungan berupa hubungan positif antar peubah x dan y serta ragam yang cukup besar. Terdapat indikasi pencilan.

4. Pembentukan Model (Fungsi lm())

Ringkasan Model Regresi

reg <- lm(y~x1+x2+x4, data = df)
summary(reg)

## 
## Call:
## lm(formula = y ~ x1 + x2 + x4, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -45.799 -16.805  -1.352   9.696  51.355 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 122.1006    11.9983  10.177 1.48e-11 ***
## x1            0.9242     0.8127   1.137   0.2639    
## x2            0.6180     0.9055   0.682   0.4999    
## x4            3.9523     1.8119   2.181   0.0366 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26.97 on 32 degrees of freedom
## Multiple R-squared:  0.386,  Adjusted R-squared:  0.3284 
## F-statistic: 6.706 on 3 and 32 DF,  p-value: 0.001223

Berdasarkan ringkasan model regresi, persamaan regresi yang terbentuk:
\[ y = 122.1+0.92x_1+0.62x_2+3.95x_4 \]
Terdapat nilai standar error residual yaitu 26.97 (kesalahan prediksinya cukup besar). Diketahui pula kualitas model berdasarkan nilai \(R^2=0.386\) yang berarti model hanya mampu menjelaskan 38.6% variasi data dan sisanya dijelaskan faktor lain di luar model. Nilai p-value = 0.0012 < 0.05, secara keseluruhan model signifikan (setidaknya ada 1 peubah bebas yang berpengaruh terhadap y). Hanya peubah x4 yang berpengaruh signifikan terhadap y, p-value = 0.037 < 0.05.

ANOVA

anova(reg)

## Analysis of Variance Table
## 
## Response: y
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## x1         1  8501.9  8501.9 11.6922 0.001731 **
## x2         1  2666.8  2666.8  3.6674 0.064461 . 
## x4         1  3459.9  3459.9  4.7582 0.036619 * 
## Residuals 32 23268.7   727.1                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan tabel ANOVA, jumlah kuadrat galat atau besarnya kesalahan total dalam model cukup besar yaitu 23268.7. Selain itu, hasil uji F menunjukkan peubah x1 dan x4 yang signifikan dengan p-value < 0.05.

Selang Kepercayaan

confint(reg)

##                  2.5 %     97.5 %
## (Intercept) 97.6609117 146.540358
## x1          -0.7310951   2.579564
## x2          -1.2265068   2.462407
## x4           0.2616205   7.643053

Hasil selang kepercayaan 95% untuk koefisien regresi menunjukkan bahwa peubah x1 dan x2 tidak terbukti berpengaruh secara signifikan terhadap y, karena rentangnya mencakup nilai nol. Peubah x4 seluruh rentangnya di atas nol atau dengan kata lain berdampak signifikan.

5. Pembentukan Model (Perhitungan Matriks Manual)

Parameter Regresi

Rumus dugaan koefisien regresi metode OLS:
\[ b_{(k+1)\times1}=(X'X)_{(k+1)\times(k+1)}^{-1}X'_{(k+1)\times{n}}y_{(n\times1)} \]

#Koef. regresi dengan metode Ordinary Least Square
x0 <- rep(1, length(y))
y <- as.matrix(y)
X <- as.matrix(cbind(x0,x1,x2,x4))
b <- solve(t(X)%*%X)%*%t(X)%*%y;round(b,4)

##        [,1]
## x0 122.1006
## x1   0.9242
## x2   0.6180
## x4   3.9523

b0 <- b[1]
b1 <- b[2]
b2 <- b[3]
b3 <- b[4]
print(paste("b0:",round(b[1],4)))

## [1] "b0: 122.1006"

print(paste("b1:",round(b[2],4)))

## [1] "b1: 0.9242"

print(paste("b2:",round(b[3],4)))

## [1] "b2: 0.618"

print(paste("b3:",round(b[4],4)))

## [1] "b3: 3.9523"

Koefisien Determinasi dan Penyesuaiannya

Rumus dugaan ragam galat:
\[ \hat{\sigma}^2=\frac{SSR_{esidual}}{n-p}=\frac{y'y-\hat{\beta}x'y}{n-p} \]

sigma_kuadrat <- (t(y)%*%y-t(b)%*%t(X)%*%y)/(n-p)
Res_se <- sqrt(sigma_kuadrat)
round(Res_se,2)

##       [,1]
## [1,] 26.97

Nilai koefisien regresi menunjukkan bahwa, secara rata-rata, prediksi model menyimpang sekitar 26.97 dari nilai aktual y.

y_duga <- b0+b1*x1+b2*x2+b3*x4
Y <- data.frame(y,y_duga)
head(Y)

##        y   y_duga
## 1 158.76 166.5139
## 2 197.19 207.5768
## 3 144.73 179.0114
## 4 140.06 173.5212
## 5 129.71 148.4538
## 6 162.59 180.6730

R_squared <- (cor(y,y_duga))^2
round(R_squared,4)

##       [,1]
## [1,] 0.386

\(R^2\) mengukur proporsi variasi dalam peubah y yang dapat dijelaskan oleh peubah bebas x dalam model. Nilai \(R^2\) sebesar 0.386 menunjukkan bahwa hanya sekitar 38.6% dari variasi dalam y dapat dijelaskan oleh model regresi. Ini menunjukkan bahwa model memiliki kecocokan yang cukup rendah dengan data.

R_squared_adj <- 1-((1-R_squared)*(n-1)/(n-p))
round(R_squared_adj,4)

##        [,1]
## [1,] 0.3284

Adjusted \(R^2\) memperhitungkan jumlah peubah bebas dalam model dan memberikan penalti untuk model yang terlalu kompleks. Nilai Adjusted \(R^2\) sebesar 0.3284 menunjukkan bahwa dengan penalti untuk jumlah peubah, model hanya mampu menjelaskan 32.84% dari variasi dalam y.

Uji F dan Standar Error Parameter Regresi

#Uji F
galat <- y-(b0+b1*x1+b2*x2+b3*x4)
KTReg <- sum((y_duga-mean(y))^2)/(p-1)

KTG <- sum(galat^2)/(n-p)

Fhit <- KTReg/KTG

dbreg <- p-1

dbg <- n-p

pv <- pf(Fhit, dbreg, dbg, lower.tail = F)

print(paste("KTR:",KTReg))

## [1] "KTR: 4876.19884277776"

print(paste("KTG",KTG))

## [1] "KTG 727.147978975695"

print(paste("Fhit:",round(Fhit,4)))

## [1] "Fhit: 6.7059"

print(paste("dbr:",dbreg))

## [1] "dbr: 3"

print(paste("dbg:",dbg))

## [1] "dbg: 32"

print(paste("pval:",pv))

## [1] "pval: 0.00122320457925806"

#Standar error koefisien regresi
se_b <- sqrt(sigma_kuadrat[1]*solve(t(X)%*%X))

## Warning in sqrt(sigma_kuadrat[1] * solve(t(X) %*% X)): NaNs produced

se_b

##           x0        x1        x2       x4
## x0 11.998293       NaN 1.1829188      NaN
## x1       NaN 0.8126577       NaN      NaN
## x2  1.182919       NaN 0.9055069      NaN
## x4       NaN       NaN       NaN 1.811898

# standar error koefisien regresi
se_b0 <- se_b[1,1];print(paste("se_b0:",round(se_b0,4)))

## [1] "se_b0: 11.9983"

se_b1 <- se_b[2,2];print(paste("se_b1:",round(se_b1,4)))

## [1] "se_b1: 0.8127"

se_b2 <- se_b[3,3];print(paste("se_b2:",round(se_b2,4)))

## [1] "se_b2: 0.9055"

se_b3 <- se_b[4,4];print(paste("se_b3:",round(se_b3,4)))

## [1] "se_b3: 1.8119"

Signifikansi Parameter (nilai t)

# t-value
t_b0 <- b0/se_b0 
t_b1 <- b1/se_b1 
t_b2 <- b2/se_b2
t_b3 <- b2/se_b3
print(paste("t_b0:",round(t_b0,2)))

## [1] "t_b0: 10.18"

print(paste("t_b1:",round(t_b1,2)))

## [1] "t_b1: 1.14"

print(paste("t_b2:",round(t_b2,2)))

## [1] "t_b2: 0.68"

print(paste("t_b3:",round(t_b3,2)))

## [1] "t_b3: 0.34"

# p-value
print(paste("p_b0:",2*pt(-abs(t_b0),db <- n-p)))

## [1] "p_b0: 1.47506344832897e-11"

print(paste("p_b1:",2*pt(-abs(t_b1),db <- n-p)))

## [1] "p_b1: 0.263858037514289"

print(paste("p_b2:",2*pt(-abs(t_b2),db <- n-p)))

## [1] "p_b2: 0.499874045483506"

print(paste("p_b3:",2*pt(-abs(t_b3),db <- n-p)))

## [1] "p_b3: 0.735294176287717"

Selang Kepercayaan (1-alpha)*100%

t <- qt(.975, db <- n-p)

BB_b0 <- b0-t*se_b0
BA_b0 <- b0+t*se_b0

BB_b1 <- b1-t*se_b1
BA_b1 <- b1+t*se_b1

BB_b2 <- b2-t*se_b2
BA_b2 <- b2+t*se_b2

BB_b3 <- b3-t*se_b3
BA_b3 <- b3+t*se_b3

Batas.Bawah <- as.matrix(c(round(BB_b0,6),round(BB_b1,6),round(BB_b2,6),round(BB_b3,6)))
Batas.Atas <- as.matrix(c(round(BA_b0,6),round(BA_b1,6),round(BA_b2,6),round(BA_b3,6)))
Fit <- (c(round(b0,6),round(b1,6),round(b2,6),round(b3,6)))

Selang.Kepercayaan <- cbind(Batas.Bawah, Fit, Batas.Atas)
colnames(Selang.Kepercayaan ) <- c("Batas Bawah Selang (2.5%)", "Fit", "Batas Atas Selang (97.5%)")
rownames(Selang.Kepercayaan ) <- c("Intersep", "b1", "b2", "b3")
Selang.Kepercayaan

##          Batas Bawah Selang (2.5%)        Fit Batas Atas Selang (97.5%)
## Intersep                 97.660912 122.100635                146.540358
## b1                       -0.731095   0.924235                  2.579564
## b2                       -1.226507   0.617950                  2.462407
## b3                        0.261620   3.952337                  7.643053

6. Perbandingan Hasil Perhitungan Manual vs Fungsi lm()

koef <- as.matrix(reg$coefficients)
penduga <- cbind(b, koef)
colnames(penduga) <- c('matriks', 'fungsi lm')
rownames(penduga) <- c("intersep", "b1", "b2", "b3")
penduga

##              matriks   fungsi lm
## intersep 122.1006350 122.1006350
## b1         0.9242345   0.9242345
## b2         0.6179503   0.6179503
## b3         3.9523366   3.9523366

Berdasar perbandingan ini diperoleh hasil yang sama yang artinya penggunaan matriks secara manual sudah tepat dalam memodelkan regresi linear berganda bergitu pula sebaliknya.

7. Perbandingan 2 Model

# Model 1
a <- summary(reg)
r.sq1 <- a$r.squared
adj_r.sq1 <- a$adj.r.squared
se1 <- a$sigma

# Model 2
df <- as.data.frame(df)
reg2 <- lm(formula = y~.-1, data = df)
b <- summary(reg2)
r.sq2 <- b$r.squared
adj_r.sq2 <- b$adj.r.squared
se2 <- b$sigma

## Membandingkan kebaikan model
model1 <- as.matrix(c(r.sq1,adj_r.sq1,se1))
model2 <- as.matrix(c(r.sq2,adj_r.sq2,se2))

tabel <- cbind(model1, model2)
colnames(tabel) <- c("Model 1", "Model 2")
rownames(tabel) <- c("R-Square", "Adj R-Square", "Standar Error Sisaan")
tabel

##                        Model 1    Model 2
## R-Square              0.386006  0.9082997
## Adj R-Square          0.328444  0.8999633
## Standar Error Sisaan 26.965682 54.6540148

Perbandingan yang dilakukan memberikan gambaran bagaimana membandingkan 2 buah model regresi berganda. Dari perbandingan ini diperoleh bahwa penggunaan Adjusted\(R^2\) akan lebih baik dari penggunaan \(R^2\) untuk membandingkan mana model terbaik.