Pada Chapter ini, kita akan membahas cara membangun model regresi menggunakan R. Langkah-langkah dapat dilakukan sebagai berikut :

  1. Pengambilan data Adapun data yang digunakan pada chapter ini adalah data tableregresisederhana.xlsx di URL https://docs.google.com/spreadsheets/d/1N5gYt9oiE42FpFryp-JGmGXvdgOD_GyP/edit?usp=sharing&ouid=117561626746546590770&rtpof=true&sd=true

Tipe file XSLX adalah merupakan tipe file dari excel 2007 . Untuk cara importnya seperti berikut;

# perintah pertama untuk membaca file data tableregresisederhana.xlsx
library(readxl)
dataregresi <- read_excel("tableregresisederhana.xlsx")

Untuk melihat dataregresi

head(dataregresi)
## # A tibble: 6 x 14
##      A1    A2    A3    A4    A5     A    B1    B2    B3     B    C1    C2    C3
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1     5     5     6     6     6    28     6     6     6    18     6     6     6
## 2     3     2     2     2     3    12     4     4     3    11     2     4     3
## 3     2     2     2     3     3    12     3     4     4    11     5     5     4
## 4     3     3     2     3     3    14     4     3     3    10     3     4     3
## 5     6     6     5     6     6    29     6     6     6    18     6     6     5
## 6     3     3     1     2     4    13     3     3     3     9     3     3     3
## # ... with 1 more variable: C <dbl>
library(ggplot2)
# Building histogram
ggplot(data=dataregresi, aes(dataregresi$A1)) +
  geom_histogram(aes(y =..density..), fill = "orange") +
  geom_density()
## Warning: Use of `dataregresi$A1` is discouraged. Use `A1` instead.

## Warning: Use of `dataregresi$A1` is discouraged. Use `A1` instead.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# loading psych package
library(psych)
## Warning: package 'psych' was built under R version 4.1.2
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
psych::describe(dataregresi)
##    vars   n  mean   sd median trimmed  mad min max range  skew kurtosis   se
## A1    1 145  4.28 1.22      4    4.31 1.48   1   6     5 -0.35    -0.36 0.10
## A2    2 145  4.23 1.22      4    4.28 1.48   1   6     5 -0.35    -0.37 0.10
## A3    3 145  4.07 1.28      4    4.14 1.48   1   6     5 -0.48    -0.32 0.11
## A4    4 145  4.08 1.31      4    4.15 1.48   1   6     5 -0.32    -0.46 0.11
## A5    5 145  4.17 1.34      4    4.27 1.48   1   6     5 -0.57    -0.17 0.11
## A     6 145 20.83 5.35     22   21.04 4.45   5  30    25 -0.54     0.31 0.44
## B1    7 145  4.47 1.05      4    4.50 1.48   1   6     5 -0.20    -0.31 0.09
## B2    8 145  4.54 1.07      5    4.57 1.48   1   6     5 -0.37    -0.44 0.09
## B3    9 145  4.53 1.14      5    4.59 1.48   1   6     5 -0.44    -0.51 0.09
## B    10 145 13.54 2.95     14   13.62 2.97   3  18    15 -0.36    -0.02 0.24
## C1   11 145  4.46 1.15      4    4.51 1.48   1   6     5 -0.43    -0.22 0.10
## C2   12 145  4.57 1.17      5    4.65 1.48   1   6     5 -0.44    -0.49 0.10
## C3   13 145  4.57 1.18      5    4.65 1.48   1   6     5 -0.53    -0.30 0.10
## C    14 145 13.60 3.25     14   13.78 2.97   3  18    15 -0.52    -0.13 0.27

Dalam contoh ini kita akan menyesuaikan model regresi menggunakan mtcars dataset R bawaan dan kemudian menghasilkan tiga plot residual yang berbeda untuk menganalisis residual.

Langkah 1: Sesuaikan model regresi.

Pertama, kita akan menyesuaikan model regresi menggunakan A sebagai variabel respon dan A1 dan A2 sebagai variabel penjelas:

#fit a regression model
model <- lm(dataregresi$A~dataregresi$A1+dataregresi$A2)

#get list of residuals 
res <- resid(model)

Langkah 2: Menghasilkan plot residual vs fitted.

Selanjutnya, kami akan menghasilkan plot residual vs fitted, yang berguna untuk mendeteksi heteroskedastisitas secara visual – mis. perubahan sistematis dalam penyebaran residu pada rentang nilai.

#produce residual vs. fitted plot
plot(fitted(model), res)

#add a horizontal line at 0 
abline(0,0)

Sumbu x menampilkan nilai fitted dan sumbu y menampilkan residu. Dari plot kita dapat melihat bahwa penyebaran residu cenderung lebih tinggi untuk nilai yang dipasang lebih tinggi, tetapi tidak terlihat cukup serius sehingga kita perlu membuat perubahan apa pun pada model.

Langkah 3: Buat plot Q-Q.

Kita juga dapat menghasilkan plot Q-Q, yang berguna untuk menentukan apakah residual mengikuti distribusi normal. Jika nilai data dalam plot jatuh sepanjang garis lurus kira-kira pada sudut 45 derajat, maka data terdistribusi normal.

#create Q-Q plot for residuals
qqnorm(res)

#add a straight diagonal line to the plot
qqline(res)

Kita dapat melihat bahwa residu cenderung menyimpang dari garis sedikit di dekat ekor, yang dapat menunjukkan bahwa mereka tidak terdistribusi secara normal.

Langkah 4: Buat plot kepadatan.

Kami juga dapat menghasilkan plot kepadatan, yang juga berguna untuk memeriksa secara visual apakah residu terdistribusi normal atau tidak. Jika plotnya kira-kira berbentuk lonceng, maka residu kemungkinan mengikuti distribusi normal.

#Create density plot of residuals
plot(density(res))

Kita dapat melihat bahwa plot kepadatan secara kasar mengikuti bentuk lonceng, meskipun sedikit miring ke kanan. Tergantung pada jenis studinya, seorang peneliti mungkin atau mungkin tidak memutuskan untuk melakukan transformasi pada data untuk memastikan bahwa residual lebih terdistribusi secara normal.

model
## 
## Call:
## lm(formula = dataregresi$A ~ dataregresi$A1 + dataregresi$A2)
## 
## Coefficients:
##    (Intercept)  dataregresi$A1  dataregresi$A2  
##          2.426           1.670           2.664
names(model)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"
plot(model$model)

model$coef
##    (Intercept) dataregresi$A1 dataregresi$A2 
##       2.426111       1.669509       2.664141
attributes(model)
## $names
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"        
## 
## $class
## [1] "lm"
round(summary(model)$coef, 3)
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)       2.426      0.702   3.458    0.001
## dataregresi$A1    1.670      0.207   8.060    0.000
## dataregresi$A2    2.664      0.208  12.817    0.000
# Pull the coefficients table from summary(lm)
crime.lm.coef <- round(summary(model)$coef, 3)
# See what this gives
class(crime.lm.coef)
## [1] "matrix" "array"
attributes(crime.lm.coef)
## $dim
## [1] 3 4
## 
## $dimnames
## $dimnames[[1]]
## [1] "(Intercept)"    "dataregresi$A1" "dataregresi$A2"
## 
## $dimnames[[2]]
## [1] "Estimate"   "Std. Error" "t value"    "Pr(>|t|)"
plot(model)

diamonds.lm <- lm(A ~ A1 + A2 + A3 + A4 + A5, data = dataregresi)

plot(diamonds.lm)

Daftar pustaka : https://www.statology.org/residual-plot-r/