1 Pendahuluan

1.1 Latar Belakang

Analisis statistika modern tidak dapat dipisahkan dari pendekatan komputasi. Melalui pemanfaatan perangkat lunak statistik seperti R, proses pemodelan, estimasi parameter, serta evaluasi model dapat dilakukan secara lebih efisien dan sistematis.

Salah satu aplikasi penting dari metode statistika adalah analisis faktor-faktor yang mempengaruhi indikator kesehatan suatu negara. Dalam konteks global, angka harapan hidup (life expectancy) merupakan indikator penting untuk menggambarkan kondisi kesehatan dan kesejahteraan suatu populasi.

Data mengenai angka harapan hidup sering digunakan untuk membandingkan tingkat pembangunan sosial dan ekonomi antar negara. Oleh karena itu, analisis statistika terhadap faktor-faktor yang mempengaruhi angka harapan hidup menjadi penting untuk memahami hubungan antar variabel yang relevan.

1.2 Tujuan Analisis

Tujuan dari analisis ini adalah:

Mengembangkan model regresi untuk menjelaskan variasi angka harapan hidup antar negara.
Melakukan evaluasi asumsi model regresi melalui analisis residual.
Mengidentifikasi kemungkinan adanya multikolinearitas antar variabel penjelas.
Mengembangkan model alternatif yang lebih baik berdasarkan evaluasi model awal.

2 Tinjauan Data

Data yang digunakan berasal dari World Health Organization (WHO) dan United Nations Population Division.

Variabel yang digunakan dalam analisis ini adalah:

LifeExpectancy2015 : angka harapan hidup saat lahir pada tahun 2015
IncomeGroup : kelompok pendapatan negara
Polio : tingkat imunisasi polio pada anak usia 1 tahun
IncomeCompositionResources : indeks komposisi pendapatan dalam Human Development Index
Schooling : rata-rata lama sekolah

Variabel LifeExpectancy2015 digunakan sebagai variabel respon, sedangkan variabel lainnya digunakan sebagai variabel penjelas.

3 Persiapan Analisis

3.1 Load Packages

Library di atas digunakan untuk:

manipulasi data (dplyr)
visualisasi (ggplot2)
analisis hubungan antar variabel (GGally)
diagnostik regresi (car)

4 Analisis Regresi

4.1 Import dan Persiapan Data

Dataset kemudian dipilih hanya pada variabel yang relevan dan observasi yang memiliki nilai hilang dihapus agar analisis dapat dilakukan secara konsisten.

4.2 Pair Plot Variabel

ggpairs(LifeExp)

Pair plot digunakan untuk mengevaluasi hubungan antar variabel serta mengidentifikasi kemungkinan pola hubungan non-linear maupun indikasi multikolinearitas.

5 Transformasi Variabel

LifeExp_transformed <- LifeExp %>%
  transmute(IncomeGroup = IncomeGroup,
            sq_Polio = Polio^2,
            IncomeCompositionResources = IncomeCompositionResources,
            Schooling = Schooling,
            LifeExpectancy2015 = LifeExpectancy2015)

ggpairs(LifeExp_transformed)

Transformasi dilakukan untuk memperbaiki hubungan linear antara variabel penjelas dan variabel respon serta untuk menstabilkan varians.

LifeExp_transformed <- LifeExp %>%
  transmute(IncomeGroup = IncomeGroup,
            exp_Polio = exp(Polio),
            IncomeCompositionResources = IncomeCompositionResources,
            Schooling = Schooling,
            LifeExpectancy2015 = LifeExpectancy2015)

ggpairs(LifeExp_transformed)

LifeExp_transformed <- LifeExp %>%
  transmute(IncomeGroup = IncomeGroup,
            Polio_6 = Polio^6,
            IncomeCompositionResources = IncomeCompositionResources,
            Schooling = Schooling,
            LifeExpectancy2015 = LifeExpectancy2015)

ggpairs(LifeExp_transformed)

Berdasarkan evaluasi grafik, transformasi Polio^6 memberikan hubungan yang relatif lebih baik dengan variabel respon dibandingkan transformasi lainnya.

6 Pembentukan Model Regresi

lm_fit <- lm(LifeExpectancy2015 ~
               IncomeGroup + Polio_6 + IncomeCompositionResources + Schooling,
             data = LifeExp_transformed)

LifeExp_transformed <- LifeExp_transformed %>%
  mutate(resid = residuals(lm_fit))

Model regresi linear digunakan untuk menjelaskan hubungan antara variabel respon dengan beberapa variabel penjelas.

7 Diagnostik Model

7.1 Plot Residual

sp1 <- ggplot(data = LifeExp_transformed,
              mapping = aes(x = resid, color = IncomeGroup))+
  geom_density()

sp2 <- ggplot(data = LifeExp_transformed,
              mapping = aes(x = Polio_6, y = resid))+
  geom_point()

sp3 <- ggplot(data = LifeExp_transformed,
              mapping = aes(x = IncomeCompositionResources, y = resid))+
  geom_point()

sp4 <- ggplot(data = LifeExp_transformed,
              mapping = aes(x = Schooling, y = resid))+
  geom_point()

sp5 <- ggplot(data = LifeExp_transformed, mapping = aes(x = resid))+
  geom_density()

grid.arrange(sp1, sp2, sp3, sp4, sp5, ncol=2)

## Warning: <ggplot> %+% x was deprecated in ggplot2 4.0.0.
## ℹ Please use <ggplot> + x instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Plot residual digunakan untuk mengevaluasi asumsi linearitas dan homoskedastisitas.

7.2 Residual vs Fitted

plot(fitted(lm_fit), LifeExp_transformed$resid)
abline(0,0, col="blue")

Plot ini membantu mengevaluasi apakah terdapat pola tertentu pada residual.

7.3 Distribusi Residual

plot(density(LifeExp_transformed$resid))
abline(v=0, col="blue")

Distribusi residual diperiksa untuk mengevaluasi asumsi normalitas.

7.4 Leverage dan Studentized Residual

car::influenceIndexPlot(lm_fit,
                        vars = c("Studentized", "hat"))

2 * length(coef(lm_fit)) / nrow(LifeExp_transformed)

## [1] 0.08235294

Observasi dengan leverage tinggi perlu diperiksa untuk memastikan tidak memberikan pengaruh berlebihan terhadap model.

8 Analisis Multikolinearitas

vif(lm_fit)

##                                 GVIF Df GVIF^(1/(2*Df))
## IncomeGroup                 6.273654  3        1.358064
## Polio_6                     1.508758  1        1.228315
## IncomeCompositionResources 14.679236  1        3.831349
## Schooling                   6.575440  1        2.564262

Nilai Variance Inflation Factor (VIF) digunakan untuk mendeteksi adanya multikolinearitas antar variabel penjelas.

Apabila nilai VIF cukup tinggi, maka variabel dengan nilai VIF terbesar dapat dipertimbangkan untuk dihapus dari model.

9 Model Alternatif

lm_new <- lm(LifeExpectancy2015 ~
               IncomeGroup + Polio_6 + Schooling,
             data = LifeExp_transformed)

Model baru dibangun dengan menghapus variabel yang memiliki indikasi multikolinearitas tinggi.

summary(lm_new)

## 
## Call:
## lm(formula = LifeExpectancy2015 ~ IncomeGroup + Polio_6 + Schooling, 
##     data = LifeExp_transformed)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.1201  -1.8831   0.4027   2.9846   8.0428 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              5.616e+01  2.871e+00  19.559  < 2e-16 ***
## IncomeGroupLow          -7.688e+00  1.596e+00  -4.818 3.28e-06 ***
## IncomeGroupLower middle -4.543e+00  1.153e+00  -3.942 0.000119 ***
## IncomeGroupUpper middle -2.827e+00  8.895e-01  -3.178 0.001773 ** 
## Polio_6                  5.199e-12  1.114e-12   4.666 6.34e-06 ***
## Schooling                1.219e+00  1.842e-01   6.621 4.84e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.012 on 164 degrees of freedom
## Multiple R-squared:  0.7523, Adjusted R-squared:  0.7447 
## F-statistic: 99.61 on 5 and 164 DF,  p-value: < 2.2e-16

Hasil estimasi menunjukkan bahwa variabel IncomeGroup, Polio_6, dan Schooling memiliki hubungan yang signifikan dengan variabel respon.

10 Deteksi Outlier

LifeExp_transformed <- LifeExp_transformed %>%
  mutate(obs_index = row_number(),
         h = hatvalues(lm_new),
         studres = rstudent(lm_new))

ggplot(data = LifeExp_transformed,
       mapping = aes(x = obs_index, y = studres))+
  geom_point() +
  geom_hline(yintercept = 0, col = "blue")

Analisis ini digunakan untuk mengidentifikasi kemungkinan adanya outlier atau observasi berpengaruh.

11 Polynomial Regression

11.1 Import Data

paper <- read_excel("Hardwood.xlsx")

Dataset ini berisi data mengenai kekuatan tarik kertas dan persentase kayu keras dalam pulp yang digunakan.

11.2 Model Polynomial Degree 2

model_d2 <- lm(tensile ~ poly(hardwood, 2, raw = TRUE),
               data = paper)

summary(model_d2)

## 
## Call:
## lm(formula = tensile ~ poly(hardwood, 2, raw = TRUE), data = paper)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.8503 -3.2482 -0.7267  4.1350  6.5506 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -6.67419    3.39971  -1.963   0.0673 .  
## poly(hardwood, 2, raw = TRUE)1 11.76401    1.00278  11.731 2.85e-09 ***
## poly(hardwood, 2, raw = TRUE)2 -0.63455    0.06179 -10.270 1.89e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.42 on 16 degrees of freedom
## Multiple R-squared:  0.9085, Adjusted R-squared:  0.8971 
## F-statistic: 79.43 on 2 and 16 DF,  p-value: 4.912e-09

ggplot(data = paper, mapping = aes(x = hardwood, y = tensile))+
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ poly(x, 2, raw = TRUE))

11.3 Model Polynomial Degree 3

model_d3 <- lm(tensile ~ poly(hardwood, 3, raw = TRUE),
               data = paper)

summary(model_d3)

## 
## Call:
## lm(formula = tensile ~ poly(hardwood, 3, raw = TRUE), data = paper)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6250 -1.6109  0.0413  1.5892  5.0216 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     5.648395   2.954663   1.912   0.0752 .  
## poly(hardwood, 3, raw = TRUE)1  3.578489   1.565854   2.285   0.0373 *  
## poly(hardwood, 3, raw = TRUE)2  0.653635   0.231330   2.826   0.0128 *  
## poly(hardwood, 3, raw = TRUE)3 -0.055188   0.009789  -5.638 4.72e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.585 on 15 degrees of freedom
## Multiple R-squared:  0.9707, Adjusted R-squared:  0.9648 
## F-statistic: 165.4 on 3 and 15 DF,  p-value: 1.025e-11

ggplot(data = paper, mapping = aes(x = hardwood, y = tensile))+
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ poly(x, 3, raw = TRUE))

Model dengan derajat 3 menunjukkan nilai RSS yang lebih kecil sehingga memberikan kecocokan model yang lebih baik.

12 Estimasi Parameter Menggunakan Matriks

X <- model.matrix(model_d2)
y <- matrix(paper$tensile)

(beta_hat <- solve( t(X) %*% X) %*% t(X) %*% y)

##                                      [,1]
## (Intercept)                    -6.6741916
## poly(hardwood, 2, raw = TRUE)1 11.7640057
## poly(hardwood, 2, raw = TRUE)2 -0.6345492

(y_hat <- X %*% beta_hat)

##         [,1]
## 1   4.455265
## 2   9.544081
## 3  14.315623
## 4  22.906883
## 5  30.229044
## 6  33.414213
## 7  36.282107
## 8  38.832727
## 9  41.066072
## 10 42.982143
## 11 44.580939
## 12 46.826707
## 13 47.803377
## 14 47.510948
## 15 45.949421
## 16 43.118796
## 17 39.019072
## 18 33.650250
## 19 27.012330

Perhitungan ini menunjukkan bahwa estimasi parameter regresi dapat diperoleh secara langsung menggunakan operasi matriks berdasarkan persamaan Ordinary Least Squares (OLS).

13 Kesimpulan

Berdasarkan hasil analisis yang telah dilakukan, model regresi yang dikembangkan mampu menjelaskan hubungan antara variabel penjelas dengan variabel respon secara memadai.

Evaluasi diagnostik menunjukkan bahwa model telah memenuhi sebagian besar asumsi regresi linear, meskipun terdapat beberapa indikasi multikolinearitas pada model awal.

Setelah dilakukan penyesuaian model dengan menghapus variabel yang memiliki multikolinearitas tinggi, model alternatif menunjukkan hasil yang lebih stabil dan dapat digunakan sebagai model akhir dalam analisis ini.

Linear and Polynomial Regression

Sri Hartati Setia Ningrum (24050123140170)

2026-03-08