Multiple linear model

Offline di Departemen Matematika

Published

October 28, 2024

Pengantar Model Regresi Linier Berganda Model regresi linier berganda digunakan untuk memprediksi nilai suatu variabel tergantung

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_n X_n + \epsilon\]

Di mana:

\(𝛽_0\) : adalah intersep atau titik potong,

\(𝛽_1,𝛽_2,\dots,𝛽_𝑝\) adalah koefisien dari setiap variabel independen, dan \(𝜖\) adalah error atau residual.

Asumsi Model Model regresi linier berganda membutuhkan beberapa asumsi untuk menghasilkan estimasi yang valid:

Linearitas: Hubungan antara variabel dependen dan independen harus linier. Independensi: Residual harus independen satu sama lain.
Homoskedastisitas: Variansi residual harus konstan di seluruh pengamatan. Normalitas: Residual harus terdistribusi normal.

Tipe-Tipe Regresi Linear di R

library(car)

Loading required package: carData

data("Salaries")
head(Salaries)

       rank discipline yrs.since.phd yrs.service  sex salary
1      Prof          B            19          18 Male 139750
2      Prof          B            20          16 Male 173200
3  AsstProf          B             4           3 Male  79750
4      Prof          B            45          39 Male 115000
5      Prof          B            40          41 Male 141500
6 AssocProf          B             6           6 Male  97000

str(Salaries)

'data.frame':   397 obs. of  6 variables:
 $ rank         : Factor w/ 3 levels "AsstProf","AssocProf",..: 3 3 1 3 3 2 3 3 3 3 ...
 $ discipline   : Factor w/ 2 levels "A","B": 2 2 2 2 2 2 2 2 2 2 ...
 $ yrs.since.phd: int  19 20 4 45 40 6 30 45 21 18 ...
 $ yrs.service  : int  18 16 3 39 41 6 23 45 20 18 ...
 $ sex          : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 1 ...
 $ salary       : int  139750 173200 79750 115000 141500 97000 175000 147765 119250 129000 ...

Model tanpa Intercept

Pada bagian ini , model dibuat tanpa menggunakan / menentukan koefisien intercept.

model_0 <- lm(salary~ 0+yrs.since.phd, data = Salaries) 
summary(model_0)


Call:
lm(formula = salary ~ 0 + yrs.since.phd, data = Salaries)

Residuals:
    Min      1Q  Median      3Q     Max 
-151053   -3962   35818   58578  102172 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
yrs.since.phd   4069.5      104.2   39.05   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 53490 on 396 degrees of freedom
Multiple R-squared:  0.7938,    Adjusted R-squared:  0.7933 
F-statistic:  1525 on 1 and 396 DF,  p-value: < 2.2e-16

Model Linear Biasa

Pada bagian ini , model dibuat hanya dengan menggunakan 1 variabel saja.

model_1var <- lm(salary~yrs.since.phd, data = Salaries) 
summary(model_1var)


Call:
lm(formula = salary ~ yrs.since.phd, data = Salaries)

Residuals:
   Min     1Q Median     3Q    Max 
-84171 -19432  -2858  16086 102383 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    91718.7     2765.8  33.162   <2e-16 ***
yrs.since.phd    985.3      107.4   9.177   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27530 on 395 degrees of freedom
Multiple R-squared:  0.1758,    Adjusted R-squared:  0.1737 
F-statistic: 84.23 on 1 and 395 DF,  p-value: < 2.2e-16

Model Linear Berganda

Pada bagian ini , model dibuat dengan minimal 2 variabel atau lebih.

model_2var <- lm(salary~yrs.since.phd + yrs.service, data = Salaries) 
summary(model_2var)


Call:
lm(formula = salary ~ yrs.since.phd + yrs.service, data = Salaries)

Residuals:
   Min     1Q Median     3Q    Max 
-79735 -19823  -2617  15149 106149 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    89912.2     2843.6  31.620  < 2e-16 ***
yrs.since.phd   1562.9      256.8   6.086 2.75e-09 ***
yrs.service     -629.1      254.5  -2.472   0.0138 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27360 on 394 degrees of freedom
Multiple R-squared:  0.1883,    Adjusted R-squared:  0.1842 
F-statistic: 45.71 on 2 and 394 DF,  p-value: < 2.2e-16

Model Linear dengan prediktor Kategorik

df <- Salaries
df$Male <- ifelse(df$sex == 'Male', 1, 0) 
head(df)

       rank discipline yrs.since.phd yrs.service  sex salary Male
1      Prof          B            19          18 Male 139750    1
2      Prof          B            20          16 Male 173200    1
3  AsstProf          B             4           3 Male  79750    1
4      Prof          B            45          39 Male 115000    1
5      Prof          B            40          41 Male 141500    1
6 AssocProf          B             6           6 Male  97000    1

model1 <- lm(salary~yrs.since.phd + yrs.service + Male, data=df)
summary(model1)


Call:
lm(formula = salary ~ yrs.since.phd + yrs.service + Male, data = df)

Residuals:
   Min     1Q Median     3Q    Max 
-79586 -19564  -3018  15071 105898 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    82875.9     4800.6  17.264  < 2e-16 ***
yrs.since.phd   1552.8      256.1   6.062 3.15e-09 ***
yrs.service     -649.8      254.0  -2.558   0.0109 *  
Male            8457.1     4656.1   1.816   0.0701 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27280 on 393 degrees of freedom
Multiple R-squared:  0.1951,    Adjusted R-squared:  0.189 
F-statistic: 31.75 on 3 and 393 DF,  p-value: < 2.2e-16

Pada bagian ini, nilai variabel sex didefinisikan sebagai nilai numerik, di mana kategori Male diubah menjadi 1 dan kategori Female diubah menjadi 0.

library(fastDummies)
dataf <- dummy_cols(Salaries,select_columns ='sex')
head(dataf)

       rank discipline yrs.since.phd yrs.service  sex salary sex_Female
1      Prof          B            19          18 Male 139750          0
2      Prof          B            20          16 Male 173200          0
3  AsstProf          B             4           3 Male  79750          0
4      Prof          B            45          39 Male 115000          0
5      Prof          B            40          41 Male 141500          0
6 AssocProf          B             6           6 Male  97000          0
  sex_Male
1        1
2        1
3        1
4        1
5        1
6        1

model2<- lm(salary~yrs.since.phd + yrs.service + sex_Male, data=dataf)
summary(model2)


Call:
lm(formula = salary ~ yrs.since.phd + yrs.service + sex_Male, 
    data = dataf)

Residuals:
   Min     1Q Median     3Q    Max 
-79586 -19564  -3018  15071 105898 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    82875.9     4800.6  17.264  < 2e-16 ***
yrs.since.phd   1552.8      256.1   6.062 3.15e-09 ***
yrs.service     -649.8      254.0  -2.558   0.0109 *  
sex_Male        8457.1     4656.1   1.816   0.0701 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27280 on 393 degrees of freedom
Multiple R-squared:  0.1951,    Adjusted R-squared:  0.189 
F-statistic: 31.75 on 3 and 393 DF,  p-value: < 2.2e-16

Pada bagian ini, nilai variabel sex didefinisikan menjadi dummy variabel , di mana kategori Male dan Female diubah menjadi kolom dengan nilai yang tepat.

Sex	Sex_Male	Sex_Female
Male	1	0
Female	0	1
Female	0	1

#tanpa buat dummy variabel
model3 <-lm(salary~yrs.since.phd + yrs.service + sex, data=Salaries)
summary(model3)


Call:
lm(formula = salary ~ yrs.since.phd + yrs.service + sex, data = Salaries)

Residuals:
   Min     1Q Median     3Q    Max 
-79586 -19564  -3018  15071 105898 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    82875.9     4800.6  17.264  < 2e-16 ***
yrs.since.phd   1552.8      256.1   6.062 3.15e-09 ***
yrs.service     -649.8      254.0  -2.558   0.0109 *  
sexMale         8457.1     4656.1   1.816   0.0701 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27280 on 393 degrees of freedom
Multiple R-squared:  0.1951,    Adjusted R-squared:  0.189 
F-statistic: 31.75 on 3 and 393 DF,  p-value: < 2.2e-16

Pada bagian ini, nilai variabel sex akan langsung digunakan secara otomatis.

#ganti base level
Salaries$sex <- relevel(Salaries$sex, ref='Male')
model4 <-lm(salary~yrs.since.phd + yrs.service + relevel(sex, ref='Male'), data=Salaries)
summary(model4)


Call:
lm(formula = salary ~ yrs.since.phd + yrs.service + relevel(sex, 
    ref = "Male"), data = Salaries)

Residuals:
   Min     1Q Median     3Q    Max 
-79586 -19564  -3018  15071 105898 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       91333.0     2941.2  31.053  < 2e-16 ***
yrs.since.phd                      1552.8      256.1   6.062 3.15e-09 ***
yrs.service                        -649.8      254.0  -2.558   0.0109 *  
relevel(sex, ref = "Male")Female  -8457.1     4656.1  -1.816   0.0701 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 27280 on 393 degrees of freedom
Multiple R-squared:  0.1951,    Adjusted R-squared:  0.189 
F-statistic: 31.75 on 3 and 393 DF,  p-value: < 2.2e-16

Regresi Linear Polinomial

datareg <- read.csv("https://raw.githubusercontent.com/aslab-math-ui/modul-prak/refs/heads/main/semuahalaman/modulprak/2024/ganjil/molin/datasets/datareg.csv")
head(datareg)

  gaji kesenangan
1    6         14
2    9         28
3   12         50
4   14         70
5   30         89
6   35         94

plot(datareg$gaji, datareg$kesenangan)

Model Linear (Orde 1)

#Model linear (orde 1)
linearModel <- lm(kesenangan ~ gaji, data=datareg)
summary(linearModel)


Call:
lm(formula = kesenangan ~ gaji, data = datareg)

Residuals:
   Min     1Q Median     3Q    Max 
-39.34 -21.99  -2.03  23.50  35.11 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  48.4531    17.3288   2.796   0.0208 *
gaji          0.2981     0.4599   0.648   0.5331  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 28.72 on 9 degrees of freedom
Multiple R-squared:  0.0446,    Adjusted R-squared:  -0.06156 
F-statistic: 0.4201 on 1 and 9 DF,  p-value: 0.5331

#Model Kuadratik (orde 2)
datareg$gajikuadrat <- datareg$gaji^2
head(datareg)

  gaji kesenangan gajikuadrat
1    6         14          36
2    9         28          81
3   12         50         144
4   14         70         196
5   30         89         900
6   35         94        1225

modelkuadratik <- lm(kesenangan~gaji + gajikuadrat, data=datareg)
modelkuadratik <- lm(kesenangan~gaji + I(gaji^2), data=datareg)
summary(modelkuadratik)


Call:
lm(formula = kesenangan ~ gaji + I(gaji^2), data = datareg)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2484 -3.7429 -0.1812  1.1464 13.6678 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -18.25364    6.18507  -2.951   0.0184 *  
gaji          6.74436    0.48551  13.891 6.98e-07 ***
I(gaji^2)    -0.10120    0.00746 -13.565 8.38e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.218 on 8 degrees of freedom
Multiple R-squared:  0.9602,    Adjusted R-squared:  0.9502 
F-statistic: 96.49 on 2 and 8 DF,  p-value: 2.51e-06

#Model Kubik (orde 3)
datareg$gajikubik <- datareg$gaji^3
modelkubik <- lm(kesenangan~gaji + gajikuadrat + gajikubik, data=datareg)
modelkubik <- lm(kesenangan~gaji + I(gaji^2) + I(gaji^3), data=datareg)
summary(modelkubik)


Call:
lm(formula = kesenangan ~ gaji + I(gaji^2) + I(gaji^3), data = datareg)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.1251 -1.2526 -0.2991  1.4780  9.4979 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.715e+01  1.076e+01  -3.454 0.010639 *  
gaji         9.610e+00  1.482e+00   6.486 0.000338 ***
I(gaji^2)   -2.022e-01  5.052e-02  -4.001 0.005182 ** 
I(gaji^3)    9.968e-04  4.949e-04   2.014 0.083855 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.289 on 7 degrees of freedom
Multiple R-squared:  0.9748,    Adjusted R-squared:  0.964 
F-statistic: 90.26 on 3 and 7 DF,  p-value: 5.856e-06

#Compare model kuadratik sama kubik
anova(modelkuadratik,modelkubik)

Analysis of Variance Table

Model 1: kesenangan ~ gaji + I(gaji^2)
Model 2: kesenangan ~ gaji + I(gaji^2) + I(gaji^3)
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1      8 309.34                              
2      7 195.84  1     113.5 4.0567 0.08385 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#plot regresi 
plot(datareg$gaji, datareg$kesenangan)
#Mengeluarkan nilai gaji
nilaigaji <- seq(0, 60, 0.1)
#Mengeluarkan nilai prediksi dari model kuadratik
predkesenangan <- predict(modelkuadratik,list(gaji=nilaigaji, gajikuadrat=nilaigaji^2))
#Mengeluarkan Plot dengan garis
lines(nilaigaji, predkesenangan, col='red')

Model Interaksi Interaksi memungkinkan kita untuk menguji apakah efek dari satu variabel independen pada \(𝑌\) dipengaruhi oleh variabel independen lain. Model interaksi dapat dituliskan sebagai:

\[𝑌=𝛽_0 +𝛽_1 𝑋_1 + 𝛽_2𝑋_2 +𝛽_3(𝑋_1×𝑋_2)+𝜖\]

Koefisien \(𝛽_3\) menunjukkan efek interaksi antara \(𝑋_1\) dan \(𝑋_2\) .

nilai_UTS <- read.csv("https://raw.githubusercontent.com/aslab-math-ui/modul-prak/refs/heads/main/semuahalaman/modulprak/2024/ganjil/molin/datasets/nilai_UTS.csv")
head(nilai_UTS)

  nilai_UTS  IQ JenisKelamin
1        93 125         Male
2        86 120       Female
3        96 115         Male
4        81 110       Female
5        92 105         Male
6        75 100       Female

str(nilai_UTS)

'data.frame':   10 obs. of  3 variables:
 $ nilai_UTS   : int  93 86 96 81 92 75 84 77 73 74
 $ IQ          : int  125 120 115 110 105 100 95 90 85 80
 $ JenisKelamin: chr  "Male" "Female" "Male" "Female" ...

GenderMale <- nilai_UTS[which(nilai_UTS$JenisKelamin=="Male"),]
GenderFemale <- nilai_UTS[which(nilai_UTS$JenisKelamin=="Female"),]

modelMale <- lm(nilai_UTS ~ IQ, data = GenderMale)
modelFemale <- lm(nilai_UTS ~ IQ, data = GenderFemale)

#Ngeluarin plot dengan garis 
plot(nilai_UTS$IQ, nilai_UTS$nilai_UTS)
lines(GenderMale$IQ, predict(modelMale), col="red",
      lty = 1 , lwd = 2)
lines(GenderFemale$IQ, predict(modelFemale), col="blue",
      lty = 2 , lwd = 2)
legend("topleft", legend=c('Male','Female'), col=c('red','blue'),
       lty = c(1,2), lwd=c(2,2))

Transformasi Variabel Jika asumsi linearitas atau homoskedastisitas tidak terpenuhi, transformasi variabel dapat membantu. Contoh transformasi yang sering digunakan meliputi:

Logaritma untuk mengatasi ketidakseimbangan skala.
Kuadrat atau akar untuk mengurangi heteroskedastisitas.

Pengujian Model Dalam mengevaluasi model, kita biasanya menggunakan beberapa metrik:

R-squared ( \(𝑅^2\) ): Mengukur proporsi variabilitas \(𝑌\) yang dapat dijelaskan oleh model.
Adjusted \(𝑅^2\) : Mengoreksi nilai \(𝑅^2\) dengan mempertimbangkan jumlah variabel independen.
F-Test: Untuk menguji apakah setidaknya satu dari variabel independen berkontribusi signifikan terhadap prediksi \(Y\).

Open Source Exploratory Data

gunakan rawgithubusercontent untuk attrive data langsung secara online.

1. Our World in Data - https://github.com/owid

Poverty, disease, hunger, climate change, war, existential risks, and inequality: the world faces many great and terrifying problems. It is these large problems that our work at Our World in Data focuses on.

2. covid - https://github.com/aatishb/covid

The purpose of this notebook is to infer the rate at which confirmed cases of COVID-19 are growing (or were growing) in various countries.

The notebook pulls data from the John Hopkins Data Repository of global Coronavirus COVID-19 cases, and then does the following things:

List total number of confirmed cases (in countries with at least 100 cases)
Attempt to fit the time series of confirmed cases in these countries to both an exponential and a logistic function
Use these curve fits to infer doubling times (i.e., time for the number of confirmed cases to double)
If the curve fit was successful, summarize doubling times for each country

3. Center for Systems Science and Engineering - https://github.com/CSSEGISandData - https://systems.jhu.edu/research/public-health/ncov

This is the data repository for the 2019 Novel Coronavirus Visual Dashboard operated by the Johns Hopkins University Center for Systems Science and Engineering (JHU CSSE). Also, Supported by ESRI Living Atlas Team and the Johns Hopkins University Applied Physics Lab (JHU APL).