Analisis Fertility Dataset (Swiss, 1888)

1. Pendahuluan

Dataset Swiss Fertility and Socioeconomic Indicators (1888) digunakan untuk memahami bagaimana faktor sosial, ekonomi, pendidikan, dan religius memengaruhi tingkat fertilitas di Swiss pada masa Transisi Demografi.

1.1 Framing Teori

Demographic Transition Theory (DTT)

DTT menjelaskan bahwa ketika masyarakat mengalami modernisasi:

1. pendidikan meningkat,
2. mortalitas menurun,
3. peran perempuan bergeser,
4. keluarga mengecilkan ukuran kelahiran.

Dalam konteks Swiss 1888:

• pendidikan menjadi engine penurunan fertilitas,
• daerah Katolik cenderung mempertahankan fertilitas tinggi (norma religius),
• sektor agraris menjaga pola keluarga besar,
• mortalitas bayi yang tinggi menciptakan replacement effect.

---

2. Data Loading dan Deskripsi Variabel

data(swiss)
head(swiss, 10)

Deskripsi singkat variabel

• Fertility : indeks fertilitas provinsi.
• Agriculture : persentase pekerja di sektor pertanian.
• Examination : persentase laki-laki lulus ujian wajib militer.
• Education : persentase laki-laki yang mengikuti pendidikan formal.
• Catholic : persentase penduduk Katolik.
• Infant.Mortality : angka kematian bayi.

---

3. Data Cleaning & Outlier Multivariat

Rumus Jarak Mahalanobis

Sebelum pemodelan, kita memeriksa kelengkapan data dan outlier multivariat.

Rumus jarak Mahalanobis:

\[ D^2 = (x - \mu)^\top \Sigma^{-1} (x - \mu) \]

sum(is.na(swiss))

## [1] 0

sum(duplicated(swiss))

## [1] 0

X <- swiss
S <- cov(X)
center <- colMeans(X)
dist_mahal <- mahalanobis(X, center, S)
cutoff <- qchisq(0.975, df = ncol(X))
outlier_idx <- which(dist_mahal > cutoff)
outlier_idx

##    La Vallee V. De Geneve 
##           19           45

insight Data tidak mengandung nilai hilang atau duplikasi. Beberapa provinsi muncul sebagai observasi multivariat yang jauh dari pusat data, hal ini wajar mengingat perbedaan sosial-historis antar provinsi. Karena ukuran sampel terbatas dan observasi berharga, data tetap dipertahankan.

---

4. Feature engineering

Distribusi Catholic cukup terpolarisasi. untuk analisis komparatif, dibuat variabel biner:

swiss$Catholic_High <- ifelse(swiss$Catholic > 50, 1, 0)
table(swiss$Catholic_High)

## 
##  0  1 
## 29 18

insight Mengkategorikan Catholic membantu menangkap perbedaan budaya antar provinsi: provinsi mayoritas Katolik cenderung mempertahankan pola keluarga besar.

---

5. Exploratory Data Analysis (EDA)

5.1 Pairwise Relationship

GGally::ggpairs(swiss)

insight

• Education memiliki korelasi negatif kuat dengan Fertility.
• Catholic terdistribusi bimodal.
• Examination dan Education berkorelasi sedang.

---

5.2 Distribusi Catholic

ggplot(swiss, aes(Catholic)) +
  geom_histogram(fill=viridis(5)[3], bins=10) +
  theme_minimal() +
  labs(title="Distribusi Variabel Catholic", x="Persentase Katolik", y="Frekuensi")

insight Distribusi sangat tidak normal, sehingga kategori Catholic_High lebih informatif.

---

5.3 Boxplot Fertility berdasarkan Catholic_High

swiss$Catholic_High <- ifelse(swiss$Catholic > 50, 1, 0)

ggplot(swiss, aes(factor(Catholic_High), Fertility, fill=factor(Catholic_High))) +
  geom_boxplot() +
  scale_fill_viridis_d() +
  theme_minimal() +
  labs(title="Fertility berdasarkan Kategori Katolik", x="Catholic High", y="Fertility")

insight Wilayah mayoritas Katolik memiliki fertilitas yang jauh lebih tinggi.

---

6. Model Building

6.1 Baseline Model

m1 <- lm(Fertility ~ Education, data=swiss)
summary(m1)

## 
## Call:
## lm(formula = Fertility ~ Education, data = swiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.036  -6.711  -1.011   9.526  19.689 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  79.6101     2.1041  37.836  < 2e-16 ***
## Education    -0.8624     0.1448  -5.954 3.66e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.446 on 45 degrees of freedom
## Multiple R-squared:  0.4406, Adjusted R-squared:  0.4282 
## F-statistic: 35.45 on 1 and 45 DF,  p-value: 3.659e-07

insight Education saja sudah menjelaskan sebagian besar variasi Fertility.

---

6.2 Full Model

m2 <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data=swiss)
summary(m2)

## 
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education + 
##     Catholic + Infant.Mortality, data = swiss)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2743  -5.2617   0.5032   4.1198  15.3213 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
## Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
## Examination      -0.25801    0.25388  -1.016  0.31546    
## Education        -0.87094    0.18303  -4.758 2.43e-05 ***
## Catholic          0.10412    0.03526   2.953  0.00519 ** 
## Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared:  0.7067, Adjusted R-squared:  0.671 
## F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

insight Education tetap menjadi prediktor paling kuat. Catholic dan Infant.Mortality juga signifikan secara positif.

6.3 Interaction Model

m3 <- lm(Fertility ~ Education * Catholic_High, data = swiss)

summary(m3)

## 
## Call:
## lm(formula = Fertility ~ Education * Catholic_High, data = swiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.476  -6.238  -1.681   6.383  16.519 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              73.4494     2.3298  31.527  < 2e-16 ***
## Education                -0.5955     0.1462  -4.074 0.000195 ***
## Catholic_High            15.7685     3.7799   4.172 0.000144 ***
## Education:Catholic_High  -0.8046     0.2896  -2.778 0.008067 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.13 on 43 degrees of freedom
## Multiple R-squared:  0.604,  Adjusted R-squared:  0.5764 
## F-statistic: 21.86 on 3 and 43 DF,  p-value: 9.405e-09

insight Model ini menambahkan interaksi antara tingkat pendidikan dan kategori wilayah Katolik.

6.4 Model Selection (Stepwise AIC)

Rumus AIC:

\[ \text{AIC} = 2k - 2 \ln(\hat{L}) \]

mod_compare <- data.frame(
  Model = c("M1", "M2", "M3"),
  Adj_R2 = c(summary(m1)$adj.r.squared, summary(m2)$adj.r.squared, summary(m3)$adj.r.squared),
  AIC = c(AIC(m1), AIC(m2), AIC(m3)),
  BIC = c(BIC(m1), BIC(m2), BIC(m3))
)
mod_compare

stepAIC(m2, direction = "both")

## Start:  AIC=190.69
## Fertility ~ Agriculture + Examination + Education + Catholic + 
##     Infant.Mortality
## 
##                    Df Sum of Sq    RSS    AIC
## - Examination       1     53.03 2158.1 189.86
## <none>                          2105.0 190.69
## - Agriculture       1    307.72 2412.8 195.10
## - Infant.Mortality  1    408.75 2513.8 197.03
## - Catholic          1    447.71 2552.8 197.75
## - Education         1   1162.56 3267.6 209.36
## 
## Step:  AIC=189.86
## Fertility ~ Agriculture + Education + Catholic + Infant.Mortality
## 
##                    Df Sum of Sq    RSS    AIC
## <none>                          2158.1 189.86
## + Examination       1     53.03 2105.0 190.69
## - Agriculture       1    264.18 2422.2 193.29
## - Infant.Mortality  1    409.81 2567.9 196.03
## - Catholic          1    956.57 3114.6 205.10
## - Education         1   2249.97 4408.0 221.43

## 
## Call:
## lm(formula = Fertility ~ Agriculture + Education + Catholic + 
##     Infant.Mortality, data = swiss)
## 
## Coefficients:
##      (Intercept)       Agriculture         Education          Catholic  
##          62.1013           -0.1546           -0.9803            0.1247  
## Infant.Mortality  
##           1.0784

insight Stepwise memberikan model yang mirip dengan model lengkap. Namun, keputusan akhir diambil juga berdasarkan pertimbangan teori: variabel yang relevan secara teoritis tetap dipertahankan meski secara heuristik stepwise menghapusnya.

---

7. Residual Diagnostics

7.1 Residual Plots

Berikut empat plot residual untuk regresi linear.

par(mfrow=c(2,2))
plot(m2)

par(mfrow=c(1,1))

insight

• Residual vs Fitted: mendekati 0 dan tanpa pola (Homoskedastisitas).
• QQ Plot: Beberapa titik ekstrem (Normal).
• Scale-Location: Fluktuasi kecil (Heteroskedastitas rendah).
• Residual vs Leverage: tidak ada titik dengan Cook’s Distance tinggi (model stabil).

7.2 Variance Inflation Factor (VIF)

\[ VIF_j = \frac{1}{1 - R_j^2} \]

car::vif(m2)

##      Agriculture      Examination        Education         Catholic 
##         2.284129         3.675420         2.774943         1.937160 
## Infant.Mortality 
##         1.107542

insight Tidak ada multikolinearitas berat (VIF < 10).

---

8. Analisis Lanjutan

8.1 Partial Effects

eff <- allEffects(m2)
plot(eff, multiline=TRUE, ci.style="bands", colors=viridis(5))

insight

• Education penurun Fertility paling kuat dan konsisten.
• Catholic menunjukan hubungan positif.
• Infant.Mortality peningkatan kematian bayi menaikkan fertilitas.

---

8.2 Bootstrap Coefficients (1000 Iterasi)

Rumus Estimator Bootstrap

\[ \hat{\beta}^{*(b)} = f(X^{*(b)}, Y^{*(b)}) \]

boot_fn <- function(data, indices){
  coef(lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data=swiss[indices, ]))
}

set.seed(123)
boot_res <- boot(swiss, boot_fn, R=1000)
boot_res

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = swiss, statistic = boot_fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##       original        bias    std. error
## t1* 66.9151817  0.2851133041 10.56779480
## t2* -0.1721140 -0.0052723381  0.06525838
## t3* -0.2580082 -0.0249958235  0.25848743
## t4* -0.8709401  0.0034267750  0.22370831
## t5*  0.1041153 -0.0008429536  0.03303239
## t6*  1.0770481  0.0215939661  0.43040004

Visualisasi Bootstrap (Distribusi Koefisien)

boot_df <- as.data.frame(boot_res$t)
colnames(boot_df) <- names(coef(m2))
boot_df <- na.omit(boot_df)

boot_melt <- melt(boot_df)

ggplot(boot_melt, aes(value, fill=variable)) +
  geom_histogram(bins=35, alpha=0.7) +
  facet_wrap(~variable, scales="free") +
  scale_fill_viridis_d() +
  theme_minimal() +
  labs(title="Distribusi Bootstrap Koefisien")

insight Bootstrap menunjukkan stabilitas koefisien tiap variabel, menegaskan kekuatan efeknya.

9. Predicted vs Actual

pred <- predict(m2)

ggplot(swiss, aes(pred, Fertility)) +
  geom_point(color="#4E598C", size=2) +
  geom_smooth(method="lm", se=FALSE, color="#1B1B3A") +
  theme_minimal() +
  labs(title="Prediksi vs Nilai Aktual", x="Prediksi", y="Aktual")

insight Garis regresi cukup dekat dengan pola aktual, model cukup akurat.

---

10. Simulasi Prediksi: Education +10

Rumus perubahan prediksi untuk kenaikan Education:

\[ \Delta \hat{Y} = \beta_{Education} \cdot \Delta Education \]

new <- swiss
new$Education <- new$Education + 10

sim <- data.frame(
  Scenario = c("Original", "Education +10"),
  Predicted = c(
    mean(predict(m2)),
    mean(predict(m2, newdata=new))
  )
)

ggplot(sim, aes(Scenario, Predicted, fill=Scenario)) +
  geom_col() +
  scale_fill_viridis_d() +
  theme_minimal() +
  labs(title="Simulasi: Dampak Kenaikan Education 10 Poin", y="Rata-rata Fertility")

insight Rata-rata Fertility turun signifikan menjelaskan pendidikan memberi dampak nyata.

---

11. Kesimpulan

Core insight

• Education adalah faktor terkuat yang menurunkan Fertility.
• Catholic dan Infant.Mortality berpengaruh positif.
• Diagnostik residual menunjukkan model valid dan stabil.
• Hasil bootstrap menegaskan kekuatan dan konsistensi koefisien.
• Simulasi menunjukkan bahwa kebijakan peningkatan pendidikan efektif menurunkan tingkat fertilitas.

Social Interpretate:

Swiss 1888 sedang bergerak dari masyarakat agraris-religius menuju modern, di mana:

• Pendidikan membuka pola pikir modern,
• Norma Katolik memperlambat transisi fertilitas,
• Modernisasi menurunkan kelahiran,
• Mortalitas bayi tinggi mendorong keluarga memperbanyak anak.