1 PENDAHULUAN

1.1 Dataset Kualitas Kesehatan

background-color:#eef4f7;border-left:6px solid #4a90e2;padding:20px;border-radius:12px;line-height:1.7;margin-bottom:20px;background-color:#eef4f7;border-left:6px solid #4a90e2;padding:20px;border-radius:12px;line-height:1.7;margin-bottom:20px;

Harapan hidup merupakan indikator penting untuk menilai kualitas kesehatan dan kesejahteraan suatu negara. Nilainya dipengaruhi oleh faktor ekonomi, pendidikan, dan tingkat kematian, serta dapat berbeda antara negara developed dan developing.

Penelitian ini menggunakan data WHO tahun 2000–2015 dengan variabel Life.expectancy, BMI, dan Adult.Mortality sebagai variabel dependen, Status sebagai faktor, serta GDP, Schooling, Income.composition.of.resources, dan HIV.AIDS sebagai kovariat.

Analisis dilakukan menggunakan tiga metode:

  • MANOVA untuk mengetahui apakah terdapat perbedaan gabungan Life.expectancy, BMI, dan Adult.Mortality berdasarkan Status.
  • ANCOVA untuk mengetahui pengaruh Status terhadap Life.expectancy setelah mengontrol GDP, Schooling, dan HIV.AIDS.
  • MANCOVA untuk mengetahui pengaruh Status terhadap gabungan variabel dependen setelah mengontrol kovariat.

2 LOAD DATASET

2.1 Load Dataset

    set.seed(123)
#load dataset 
data <- read.csv("D:/coding/coding sem 4/Analisis Multivariat/tugas harian/tugas 9/manova 3.csv")
data$X <- NULL
data$IPM_Category <- cut(
  data$Income.composition.of.resources,
  breaks = c(-Inf, 0.55, 0.70, 0.80, 1),
  labels = c("Rendah", "Sedang", "Tinggi", "Sangat Tinggi"),
  right = FALSE
)

data$IPM_Category <- as.character(data$IPM_Category)
# load library
library(MVN)
library(biotools)
library(car)
head(data)
##       Country Year     Status Life.expectancy Adult.Mortality infant.deaths
## 1 Afghanistan 2015 Developing            65.0             263            62
## 2 Afghanistan 2014 Developing            59.9             271            64
## 3 Afghanistan 2013 Developing            59.9             268            66
## 4 Afghanistan 2012 Developing            59.5             272            69
## 5 Afghanistan 2011 Developing            59.2             275            71
## 6 Afghanistan 2010 Developing            58.8             279            74
##   Alcohol percentage.expenditure Hepatitis.B Measles  BMI under.five.deaths
## 1    0.01              71.279624          65    1154 19.1                83
## 2    0.01              73.523582          62     492 18.6                86
## 3    0.01              73.219243          64     430 18.1                89
## 4    0.01              78.184215          67    2787 17.6                93
## 5    0.01               7.097109          68    3013 17.2                97
## 6    0.01              79.679367          66    1989 16.7               102
##   Polio Total.expenditure Diphtheria HIV.AIDS       GDP Population
## 1     6              8.16         65      0.1 584.25921   33736494
## 2    58              8.18         62      0.1 612.69651     327582
## 3    62              8.13         64      0.1 631.74498   31731688
## 4    67              8.52         67      0.1 669.95900    3696958
## 5    68              7.87         68      0.1  63.53723    2978599
## 6    66              9.20         66      0.1 553.32894    2883167
##   thinness..1.19.years thinness.5.9.years Income.composition.of.resources
## 1                 17.2               17.3                           0.479
## 2                 17.5               17.5                           0.476
## 3                 17.7               17.7                           0.470
## 4                 17.9               18.0                           0.463
## 5                 18.2               18.2                           0.454
## 6                 18.4               18.4                           0.448
##   Schooling IPM_Category
## 1      10.1       Rendah
## 2      10.0       Rendah
## 3       9.9       Rendah
## 4       9.8       Rendah
## 5       9.5       Rendah
## 6       9.2       Rendah

2.2 keperluan Transform

#Y <- sqrt(y)
YY <- as.data.frame(scale(
  data[, c("Life.expectancy", "BMI", "Adult.Mortality")]
))

library(bestNormalize)

YYY <- data.frame(
  Life.expectancy = bestNormalize::yeojohnson(data$Life.expectancy)$x.t,
  BMI             = bestNormalize::yeojohnson(data$BMI)$x.t,
  Adult.Mortality = bestNormalize::yeojohnson(data$Adult.Mortality)$x.t
)
##Mengganti outlier dengan batas maksimum/minimum IQR
#cap_outlier_group <- function(x) {
#  Q1 <- quantile(x, 0.25, na.rm = TRUE)
#  Q3 <- quantile(x, 0.75, na.rm = TRUE)
#  IQRv <- Q3 - Q1
#
#  lower <- Q1 - 1.5 * IQRv
#  upper <- Q3 + 1.5 * IQRv

#  x[x < lower] <- lower
#  x[x > upper] <- upper
#  x
#}

2.3 Dependen(y),Faktor,dan Kovariat

# Dependen (Y)
#y <- data[, c("Life.expectancy", "BMI", "Adult.Mortality")]
y<-YYY

# Faktor
data$Status <- as.factor(data$Status)
data$IPM_Category <- as.factor(data$IPM_Category)

# Faktor yang digunakan pada two-way
factors <- data[, c("Status", "IPM_Category")]

# Kovariat
covars <- data[, c("GDP", "Schooling", "HIV.AIDS")]

data$Group <- interaction(data$Status, data$IPM_Category)

3 VISUALISASI DISTRIBUSI

Bagian ini menampilkan eksplorasi data melalui statistik deskriptif, histogram, boxplot, pie chart outlier, dan korelasi antar variabel.

3.1 Statistik

library(dplyr)
library(tidyr)

stats_table <- data %>%
  summarise(
    across(
      c(Life.expectancy, BMI, Adult.Mortality, GDP, Schooling, HIV.AIDS),
      list(
        Min    = ~min(., na.rm = TRUE),
        Q1     = ~quantile(., 0.25, na.rm = TRUE),
        Mean   = ~mean(., na.rm = TRUE),
        Median = ~median(., na.rm = TRUE),
        Q3     = ~quantile(., 0.75, na.rm = TRUE),
        Max    = ~max(., na.rm = TRUE)
      )
    )
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = c("Variable", ".value"),
    names_sep = "_"
  )

stats_table
## # A tibble: 6 × 7
##   Variable          Min    Q1    Mean Median     Q3      Max
##   <chr>           <dbl> <dbl>   <dbl>  <dbl>  <dbl>    <dbl>
## 1 Life.expectancy 36.3   63.2   69.2    72.1   75.7     89  
## 2 BMI              1     19.4   38.3    43.2   56.1     87.3
## 3 Adult.Mortality  1     74    165.    144    227      723  
## 4 GDP              1.68 459.  6745.   1633.  5368.  119173. 
## 5 Schooling        0     10     11.9    12.3   14.2     20.7
## 6 HIV.AIDS         0.1    0.1    1.74    0.1    0.8     50.6

3.2 Histogram Density

library(ggplot2)
library(gridExtra)

p1 <- ggplot(data, aes(x = Life.expectancy)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 25,
    fill = "skyblue",
    color = "black"
  ) +
  geom_density(color = "blue", linewidth = 1) +
  ggtitle("Distribusi Life Expectancy") +
  theme_minimal()

p2 <- ggplot(data, aes(x = BMI)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 25,
    fill = "lightgreen",
    color = "black"
  ) +
  geom_density(color = "darkgreen", linewidth = 1) +
  ggtitle("Distribusi BMI") +
  theme_minimal()

p3 <- ggplot(data, aes(x = Adult.Mortality)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 25,
    fill = "salmon",
    color = "black"
  ) +
  geom_density(color = "red", linewidth = 1) +
  ggtitle("Distribusi Adult Mortality") +
  theme_minimal()

grid.arrange(p1, p2, p3, ncol = 1)

3.3 Histogram Density 2

library(ggplot2)
library(gridExtra)

p1 <- ggplot(data, aes(x = Life.expectancy)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 25,
    fill = "skyblue",
    color = "black"
  ) +
  geom_density(color = "blue", linewidth = 1) +
  ggtitle("Distribusi Life Expectancy") +
  theme_minimal()

p2 <- ggplot(data, aes(x = BMI)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 25,
    fill = "lightgreen",
    color = "black"
  ) +
  geom_density(color = "darkgreen", linewidth = 1) +
  ggtitle("Distribusi BMI") +
  theme_minimal()

p3 <- ggplot(data, aes(x = Adult.Mortality)) +
  geom_histogram(
    aes(y = after_stat(density)),
    bins = 25,
    fill = "salmon",
    color = "black"
  ) +
  geom_density(color = "red", linewidth = 1) +
  ggtitle("Distribusi Adult Mortality") +
  theme_minimal()

grid.arrange(p1, p2, p3, ncol = 3)

3.4 Histogram Count

library(ggplot2)
library(gridExtra)

p1 <- ggplot(data, aes(x = Life.expectancy)) +
  geom_histogram(
    bins = 25,
    fill = "skyblue",
    color = "gray40"
  ) +
  labs(
    title = "Life Expectancy",
    x = NULL,
    y = "Count"
  ) +
  theme_minimal()

p2 <- ggplot(data, aes(x = BMI)) +
  geom_histogram(
    bins = 25,
    fill = "skyblue",
    color = "gray40"
  ) +
  labs(
    title = "BMI",
    x = NULL,
    y = NULL
  ) +
  theme_minimal()

p3 <- ggplot(data, aes(x = Adult.Mortality)) +
  geom_histogram(
    bins = 25,
    fill = "skyblue",
    color = "gray40"
  ) +
  labs(
    title = "Adult Mortality",
    x = NULL,
    y = NULL
  ) +
  theme_minimal()

grid.arrange(
  p1, p2, p3,
  ncol = 3,
  widths = c(1.3, 1.3, 1.3)
)

3.5 Boxplot

library(ggplot2)
library(gridExtra)

# Boxplot Life Expectancy
p1 <- ggplot(data, aes(x = Status, y = Life.expectancy)) +
  geom_boxplot(fill = "skyblue", alpha = 0.8) +
  labs(
    title = "(a) Life Expectancy",
    x = "Status",
    y = "Life Expectancy"
  ) +
  theme_minimal(base_size = 13)

# Boxplot BMI
p2 <- ggplot(data, aes(x = Status, y = BMI)) +
  geom_boxplot(fill = "skyblue", alpha = 0.8) +
  labs(
    title = "(b) BMI",
    x = "Status",
    y = "BMI"
  ) +
  theme_minimal(base_size = 13)

# Boxplot Adult Mortality
p3 <- ggplot(data, aes(x = Status, y = Adult.Mortality)) +
  geom_boxplot(fill = "skyblue", alpha = 0.8) +
  labs(
    title = "(c) Adult Mortality",
    x = "Status",
    y = "Adult Mortality"
  ) +
  theme_minimal(base_size = 13)

grid.arrange(p1, p2, p3, ncol = 3)

library(ggplot2)
library(gridExtra)

# Boxplot Life Expectancy
p1 <- ggplot(data, aes(x = IPM_Category, y = Life.expectancy)) +
  geom_boxplot(fill = "skyblue", alpha = 0.8) +
  labs(
    title = "(a) Life Expectancy",
    x = "IPM_Category",
    y = "Life Expectancy"
  ) +
  theme_minimal(base_size = 13,) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Boxplot BMI
p2 <- ggplot(data, aes(x = IPM_Category, y = BMI)) +
  geom_boxplot(fill = "skyblue", alpha = 0.8) +
  labs(
    title = "(b) BMI",
    x = "IPM_Category",
    y = "BMI"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

# Boxplot Adult Mortality
p3 <- ggplot(data, aes(x = IPM_Category, y = Adult.Mortality)) +
  geom_boxplot(fill = "skyblue", alpha = 0.8) +
  labs(
    title = "(c) Adult Mortality",
    x = "IPM_Category",
    y = "Adult Mortality"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

grid.arrange(p1, p2, p3, ncol = 3)

library(ggplot2)
library(gridExtra)

p1 <- ggplot(data, aes(x = Status, y = Life.expectancy, fill = IPM_Category)) +
  geom_boxplot(
    alpha = 0.8,
    position = position_dodge(width = 0.8),
    outlier.size = 1
  ) +
  labs(
    title = "(a) Life Expectancy",
    x = "Status",
    y = "Life Expectancy",
    fill = "IPM"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

p2 <- ggplot(data, aes(x = Status, y = BMI, fill = IPM_Category)) +
  geom_boxplot(
    alpha = 0.8,
    position = position_dodge(width = 0.8),
    outlier.size = 1
  ) +
  labs(
    title = "(b) BMI",
    x = "Status",
    y = "BMI",
    fill = "IPM"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

p3 <- ggplot(data, aes(x = Status, y = Adult.Mortality, fill = IPM_Category)) +
  geom_boxplot(
    alpha = 0.8,
    position = position_dodge(width = 0.8),
    outlier.size = 1
  ) +
  labs(
    title = "(c) Adult Mortality",
    x = "Status",
    y = "Adult Mortality",
    fill = "IPM"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

grid.arrange(p1, p2, p3, ncol = 3)

3.6 Pie Outlier

# Hitung outlier multivariat
YYY <- na.omit(data[, c("Life.expectancy", "BMI", "Adult.Mortality")])

md <- mahalanobis(YYY, colMeans(YYY), cov(YYY))
cutoff <- qchisq(0.975, df = 3)

outlier_idx <- which(md > cutoff)

# Jumlah outlier dan non-outlier
n_outlier <- length(outlier_idx)
n_normal  <- nrow(YYY) - n_outlier

# Data pie chart
pie_data <- data.frame(
  kategori = c("Outlier", "Non-Outlier"),
  jumlah   = c(n_outlier, n_normal)
)

library(ggplot2)

ggplot(pie_data, aes(x = "", y = jumlah, fill = kategori)) +
  geom_col(width = 1, color = "white") +
  coord_polar(theta = "y") +
  geom_text(
    aes(label = paste0(kategori, "\n", round(100 * jumlah / sum(jumlah), 1), "%")),
    position = position_stack(vjust = 0.5),
    color = "black",
    size = 5
  ) +
  labs(title = "Proporsi Outlier Multivariat") +
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

3.7 Korelasi

library(corrplot)

# Ambil hanya variabel numerik yang dipakai
vars <- data[, c(
  "Life.expectancy",
  "BMI",
  "Adult.Mortality",
  "GDP",
  "Schooling",
  "Income.composition.of.resources",
  "HIV.AIDS"
)]

# Matriks korelasi
cor_mat <- cor(vars, use = "complete.obs")

# Heatmap korelasi
corrplot(
  cor_mat,
  method = "color",
  type = "lower",
  addCoef.col = "black",
  tl.col = "black",
  tl.srt = 45,
  number.cex = 0.7
)

library(corrplot)

# Ambil hanya variabel numerik yang dipakai
library(corrplot)

# Korelasi antar variabel dependen hasil transformasi
cor_mat_y <- cor(YYY, use = "complete.obs")

corrplot(
  cor_mat_y,
  method = "color",
  type = "lower",
  addCoef.col = "black",
  tl.col = "black",
  tl.srt = 45,
  number.cex = 0.9
)

4 ASUSMSI (pakai data transform)

4.1 Korelasi

cor(y, use = "complete.obs")
##                 Life.expectancy        BMI Adult.Mortality
## Life.expectancy       1.0000000  0.5489256      -0.5881842
## BMI                   0.5489256  1.0000000      -0.3413626
## Adult.Mortality      -0.5881842 -0.3413626       1.0000000

4.2 Uji Normalitas

Uji Normalitas Multivariat dan Univariat

Henze–Zirkler Test / Uji MVN HZ

Digunakan untuk mengecek apakah data berdistribusi normal multivariat.

H0: data mengikuti distribusi normal multivariat.
H1: data tidak mengikuti distribusi normal multivariat.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → data normal multivariat
  • p-value ≤ 0,05 → tolak H0 → data tidak normal multivariat

Anderson–Darling Test / Uji Anderson–Darling

Digunakan untuk menguji apakah suatu variabel mengikuti distribusi tertentu, paling sering distribusi normal.

H0: data berdistribusi normal.
H1: data tidak berdistribusi normal.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → data berdistribusi normal
  • p-value ≤ 0,05 → tolak H0 → data tidak berdistribusi normal 
library(MVN)
result <- mvn(
  data = y,
  mvn_test = "hz",
  univariate_test = "AD",
  tidy = FALSE
)

result$multivariate_normality
##            Test Statistic p.value     Method
## 1 Henze-Zirkler  96.44042       0 asymptotic
result$univariate_normality
##               Test        Variable Statistic      p.value
## 1 Anderson-Darling Life.expectancy 18.691659 3.700000e-24
## 2 Anderson-Darling             BMI 79.251711 3.700000e-24
## 3 Anderson-Darling Adult.Mortality  9.374759 1.053285e-22

4.3 Uji Dependensi

Uji Independensi (Bartlett Test)

Bartlett’s Test of Sphericity / Uji Bartlett

Digunakan untuk melihat apakah terdapat korelasi antar variabel dependen.

H0: matriks korelasi = matriks identitas (tidak ada korelasi).
H1: terdapat korelasi antar variabel.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → tidak ada korelasi
  • p-value ≤ 0,05 → tolak H0 → ada korelasi 
library(psych)

# One-way : berdasarkan Status
cat(">>> One-Way: Status >>>>\n\n")
## >>> One-Way: Status >>>>
for(g in levels(data$Status)) {
  idx <- data$Status == g
  
  cat("Group:", g, "\n")
  print(
    cortest.bartlett(
      cor(y[idx, ], use = "complete.obs"),
      n = sum(complete.cases(y[idx, ]))
    )
  )
  cat("\n")
}
## Group: Developed 
## $chisq
## [1] 61.11845
## 
## $p.value
## [1] 3.390489e-13
## 
## $df
## [1] 3
## 
## 
## Group: Developing 
## $chisq
## [1] 1698.701
## 
## $p.value
## [1] 0
## 
## $df
## [1] 3
# One-way : berdasarkan IPM_Category
cat(">>> One-Way: IPM_Category >>>n\n")
## >>> One-Way: IPM_Category >>>n
for(g in levels(data$IPM_Category)) {
  idx <- data$IPM_Category == g
  
  cat("Group:", g, "\n")
  print(
    cortest.bartlett(
      cor(y[idx, ], use = "complete.obs"),
      n = sum(complete.cases(y[idx, ]))
    )
  )
  cat("\n")
}
## Group: Rendah 
## $chisq
## [1] 263.0804
## 
## $p.value
## [1] 9.69287e-57
## 
## $df
## [1] 3
## 
## 
## Group: Sangat Tinggi 
## $chisq
## [1] 30.6281
## 
## $p.value
## [1] 1.017973e-06
## 
## $df
## [1] 3
## 
## 
## Group: Sedang 
## $chisq
## [1] 221.4836
## 
## $p.value
## [1] 9.594094e-48
## 
## $df
## [1] 3
## 
## 
## Group: Tinggi 
## $chisq
## [1] 125.3309
## 
## $p.value
## [1] 5.484342e-27
## 
## $df
## [1] 3
# Two-way : Status × IPM_Category
cat(">>>>Two-Way: Status × IPM_Category >>>>n\n")
## >>>>Two-Way: Status × IPM_Category >>>>n
for(g in levels(data$Group)) {
  idx <- data$Group == g
  
  cat("Group:", g, "\n")
  print(
    cortest.bartlett(
      cor(y[idx, ], use = "complete.obs"),
      n = sum(complete.cases(y[idx, ]))
    )
  )
  cat("\n")
}
## Group: Developed.Rendah 
## $chisq
## [1] 1.231886
## 
## $p.value
## [1] 0.7453674
## 
## $df
## [1] 3
## 
## 
## Group: Developing.Rendah 
## $chisq
## [1] 240.8447
## 
## $p.value
## [1] 6.249352e-52
## 
## $df
## [1] 3
## 
## 
## Group: Developed.Sangat Tinggi 
## $chisq
## [1] 35.5575
## 
## $p.value
## [1] 9.288143e-08
## 
## $df
## [1] 3
## 
## 
## Group: Developing.Sangat Tinggi 
## $chisq
## [1] 1.602302
## 
## $p.value
## [1] 0.6588681
## 
## $df
## [1] 3
## 
## 
## Group: Developed.Sedang 
## $chisq
## [1] 2.300024
## 
## $p.value
## [1] 0.5125163
## 
## $df
## [1] 3
## 
## 
## Group: Developing.Sedang 
## $chisq
## [1] 188.7532
## 
## $p.value
## [1] 1.134831e-40
## 
## $df
## [1] 3
## 
## 
## Group: Developed.Tinggi 
## $chisq
## [1] 13.12114
## 
## $p.value
## [1] 0.004381821
## 
## $df
## [1] 3
## 
## 
## Group: Developing.Tinggi 
## $chisq
## [1] 109.5443
## 
## $p.value
## [1] 1.375354e-23
## 
## $df
## [1] 3
library(psych)
cor_mat <- cor(y, use = "complete.obs") 
cortest.bartlett( cor_mat, n = nrow(na.omit(y)) )
## $chisq
## [1] 2300.871
## 
## $p.value
## [1] 0
## 
## $df
## [1] 3
library(car)

model_cov <- lm(
  GDP ~ Schooling + HIV.AIDS,
  data = data
)

vif(model_cov)
## Schooling  HIV.AIDS 
##  1.052744  1.052744

4.4 Uji Homogenitas

Uji Homogenitas / Identik Univariat (Levene Test)

Levene’s Test / Uji Levene

Digunakan untuk menguji homogenitas varians antar grup.

H0: varians semua grup sama (homogen).
H1: ada minimal satu grup dengan varians berbeda.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → varians homogen
  • p-value ≤ 0,05 → tolak H0 → varians tidak homogen 
library(car)
cat("\nLife.expectancy (transformasi)\n")
## 
## Life.expectancy (transformasi)
leveneTest(YYY$Life.expectancy ~ data$Group)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    7  92.363 < 2.2e-16 ***
##       2930                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nBMI (transformasi)\n")
## 
## BMI (transformasi)
leveneTest(YYY$BMI ~ data$Group)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    7  13.695 < 2.2e-16 ***
##       2930                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nAdult.Mortality (transformasi)\n")
## 
## Adult.Mortality (transformasi)
leveneTest(YYY$Adult.Mortality ~ data$Group)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    7   90.61 < 2.2e-16 ***
##       2930                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nLife Expectancy\n")
## 
## Life Expectancy
print(
  bartlett.test(
    YYY$Life.expectancy ~ data$Group
  )
)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  YYY$Life.expectancy by data$Group
## Bartlett's K-squared = 798.23, df = 7, p-value < 2.2e-16
cat("\nBMI\n")
## 
## BMI
print(
  bartlett.test(
    YYY$BMI ~ data$Group
  )
)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  YYY$BMI by data$Group
## Bartlett's K-squared = 146.79, df = 7, p-value < 2.2e-16
cat("\nAdult Mortality\n")
## 
## Adult Mortality
print(
  bartlett.test(
    YYY$Adult.Mortality ~ data$Group
  )
)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  YYY$Adult.Mortality by data$Group
## Bartlett's K-squared = 1122.8, df = 7, p-value < 2.2e-16

Uji Homogenitas / Identik Multivariat (Box’s M Test)

Box’s M Test / Uji Box’s M

Digunakan untuk menguji kesamaan matriks kovarians antar grup pada MANOVA/MANCOVA.

H0: matriks kovarians semua grup sama.
H1: ada matriks kovarians grup yang berbeda.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → matriks kovarians homogen
  • p-value ≤ 0,05 → tolak H0 → matriks kovarians tidak homogen 
library(biotools)

boxM(y, data$Group)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  y
## Chi-Sq (approx.) = 773.07, df = 42, p-value < 2.2e-16
library(biotools)

# One-Way MANOVA : Status
cat("\nstatus\n")
## 
## status
box_status <- boxM(y, data$Status)
box_status
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  y
## Chi-Sq (approx.) = 348.08, df = 6, p-value < 2.2e-16
# One-Way MANOVA : IPM_Category
cat("\nipm\n")
## 
## ipm
box_ipm <- boxM(y, data$IPM_Category)
box_ipm
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  y
## Chi-Sq (approx.) = 833.14, df = 18, p-value < 2.2e-16
# Two-Way MANOVA : Status × IPM_Category
cat("\nstatus X ipm\n")
## 
## status X ipm
box_group <- boxM(y, data$Group)
box_group
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  y
## Chi-Sq (approx.) = 773.07, df = 42, p-value < 2.2e-16

5 MODEL

5.1 Manova 2 way

MANOVA (Multivariate Analysis of Variance)

Digunakan untuk menguji perbedaan rata-rata beberapa variabel dependen sekaligus antar grup.

H0: vektor rata-rata semua grup sama.
H1: minimal ada satu grup dengan vektor rata-rata berbeda.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → tidak ada perbedaan multivariat
  • p-value ≤ 0,05 → tolak H0 → ada perbedaan multivariat 
manova_model <- manova(
  cbind(Life.expectancy, BMI, Adult.Mortality) ~
    Status * IPM_Category,
  data = data
)

summary(manova_model, test = "Pillai")
##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## Status                 1 0.44553   784.25      3   2928 < 2.2e-16 ***
## IPM_Category           3 0.62732   258.22      9   8790 < 2.2e-16 ***
## Status:IPM_Category    3 0.03993    13.18      9   8790 < 2.2e-16 ***
## Residuals           2930                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.2 Anova

ANOVA (Analysis of Variance)

Digunakan untuk menguji apakah rata-rata satu variabel dependen berbeda antar grup.

H0: semua rata-rata grup sama.
H1: minimal ada satu rata-rata grup yang berbeda.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → tidak ada perbedaan rata-rata
  • p-value ≤ 0,05 → tolak H0 → ada perbedaan rata-rata 
summary.aov(manova_model)
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 2161.874 < 2.2e-16 ***
## IPM_Category           3 118282   39427 1385.602 < 2.2e-16 ***
## Status:IPM_Category    3   2592     864   30.363 < 2.2e-16 ***
## Residuals           2930  83373      28                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq F value    Pr(>F)    
## Status                 1 113688  113688 441.873 < 2.2e-16 ***
## IPM_Category           3 294160   98053 381.104 < 2.2e-16 ***
## Status:IPM_Category    3  10586    3529  13.715 6.985e-09 ***
## Residuals           2930 753851     257                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 419.5214 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 287.5568 < 2.2e-16 ***
## Status:IPM_Category    3   278144   92715   8.6813 9.849e-06 ***
## Residuals           2930 31291795   10680                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.3 Ancova

ANCOVA (Analysis of Covariance)

Digunakan untuk membandingkan rata-rata grup sambil mengontrol pengaruh kovariat.

H0: setelah dikontrol kovariat, rata-rata grup sama.
H1: setelah dikontrol kovariat, rata-rata grup berbeda.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → tidak ada pengaruh grup
  • p-value ≤ 0,05 → tolak H0 → ada pengaruh grup setelah kontrol kovariat 
ancova_model <- lm(
  Life.expectancy ~
    Status * IPM_Category +
    GDP + Schooling + HIV.AIDS,
  data = data
)

anova(ancova_model)
## Analysis of Variance Table
## 
## Response: Life.expectancy
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.4 Mancova

MANCOVA (Multivariate Analysis of Covariance)

Digunakan untuk menguji perbedaan beberapa variabel dependen antar grup dengan mengontrol kovariat.

H0: setelah dikontrol kovariat, vektor rata-rata semua grup sama.
H1: setelah dikontrol kovariat, minimal ada satu grup berbeda.

Keputusan:

  • p-value > 0,05 → gagal tolak H0 → tidak ada perbedaan multivariat
  • p-value ≤ 0,05 → tolak H0 → ada perbedaan multivariat setelah kontrol kovariat 
mancova_model <- manova(
  cbind(Life.expectancy, BMI, Adult.Mortality) ~
    Status * IPM_Category +
    GDP + Schooling + HIV.AIDS,
  data = data
)

summary(mancova_model, test = "Pillai")
##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## Status                 1 0.54928  1188.18      3   2925 < 2.2e-16 ***
## IPM_Category           3 0.71833   307.17      9   8781 < 2.2e-16 ***
## GDP                    1 0.00770     7.57      3   2925 4.843e-05 ***
## Schooling              1 0.06331    65.90      3   2925 < 2.2e-16 ***
## HIV.AIDS               1 0.36741   566.29      3   2925 < 2.2e-16 ***
## Status:IPM_Category    3 0.04490    14.83      9   8781 < 2.2e-16 ***
## Residuals           2927                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6 PERBANDINGAN DATA ASLI DAN DATA DATA TRANFORM

background-color:#eef4f7; border-left:6px solid #4a90e2; padding:20px; border-radius:12px; line-height:1.7; margin-bottom:20px; “} Karena data asli terlihat skew, maka dilakukan perbandingan antara data asli dan data yang telah ditransformasi menggunakan metode Yeo-Johnson. Transformasi ini bertujuan memperbaiki distribusi data agar lebih mendekati normal.

6.1 VISUALISASI

6.1.1 QQplot

library(ggplot2)
library(gridExtra)
library(bestNormalize)

# Data sebelum transformasi
Y_before <- data[, c("Life.expectancy", "BMI", "Adult.Mortality")]

# Data sesudah transformasi
Y_after <- data.frame(
  Life.expectancy = bestNormalize::yeojohnson(data$Life.expectancy)$x.t,
  BMI             = bestNormalize::yeojohnson(data$BMI)$x.t,
  Adult.Mortality = bestNormalize::yeojohnson(data$Adult.Mortality)$x.t
)

# QQ plot sebelum
p1 <- ggplot(Y_before, aes(sample = Life.expectancy)) +
  stat_qq() + stat_qq_line(color = "red") +
  ggtitle("Sebelum: Life Expectancy") +
  theme_minimal()

p2 <- ggplot(Y_before, aes(sample = BMI)) +
  stat_qq() + stat_qq_line(color = "red") +
  ggtitle("Sebelum: BMI") +
  theme_minimal()

p3 <- ggplot(Y_before, aes(sample = Adult.Mortality)) +
  stat_qq() + stat_qq_line(color = "red") +
  ggtitle("Sebelum: Adult Mortality") +
  theme_minimal()

# QQ plot sesudah
p4 <- ggplot(Y_after, aes(sample = Life.expectancy)) +
  stat_qq() + stat_qq_line(color = "blue") +
  ggtitle("Sesudah: Life Expectancy") +
  theme_minimal()

p5 <- ggplot(Y_after, aes(sample = BMI)) +
  stat_qq() + stat_qq_line(color = "blue") +
  ggtitle("Sesudah: BMI") +
  theme_minimal()

p6 <- ggplot(Y_after, aes(sample = Adult.Mortality)) +
  stat_qq() + stat_qq_line(color = "blue") +
  ggtitle("Sesudah: Adult Mortality") +
  theme_minimal()

grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3)

6.1.2 HISTogram

library(ggplot2)
library(gridExtra)
library(bestNormalize)

# Data asli
Y_before <- data[, c("Life.expectancy", "BMI", "Adult.Mortality")]

# Data transformasi
Y_after <- data.frame(
  Life.expectancy = bestNormalize::yeojohnson(data$Life.expectancy)$x.t,
  BMI             = bestNormalize::yeojohnson(data$BMI)$x.t,
  Adult.Mortality = bestNormalize::yeojohnson(data$Adult.Mortality)$x.t
)

# Fungsi histogram + density
plot_hist <- function(df, var, judul, warna) {
  ggplot(df, aes_string(x = var)) +
    geom_histogram(
      aes(y = after_stat(density)),
      bins = 25,
      fill = warna,
      color = "black",
      alpha = 0.7
    ) +
    geom_density(color = "red", linewidth = 1) +
    theme_minimal() +
    labs(title = judul, x = var, y = "Density")
}

# Sebelum
p1 <- plot_hist(Y_before, "Life.expectancy", "Sebelum: Life Expectancy", "skyblue")
p2 <- plot_hist(Y_before, "BMI", "Sebelum: BMI", "lightgreen")
p3 <- plot_hist(Y_before, "Adult.Mortality", "Sebelum: Adult Mortality", "salmon")

# Sesudah
p4 <- plot_hist(Y_after, "Life.expectancy", "Sesudah: Life Expectancy", "skyblue")
p5 <- plot_hist(Y_after, "BMI", "Sesudah: BMI", "lightgreen")
p6 <- plot_hist(Y_after, "Adult.Mortality", "Sesudah: Adult Mortality", "salmon")

grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3)

6.1.3 boxplot

library(ggplot2)
library(gridExtra)

# Tambahkan hasil transformasi YYY ke data
data_plot <- data
data_plot$Life_t <- YYY$Life.expectancy
data_plot$BMI_t  <- YYY$BMI
data_plot$Mort_t <- YYY$Adult.Mortality

# Sebelum transformasi
p1 <- ggplot(data_plot, aes(x = Status, y = Life.expectancy, fill = IPM_Category)) +
  geom_boxplot(alpha = 0.8, position = position_dodge(width = 0.8), outlier.size = 0.8) +
  labs(title = "(a) Life Expectancy - Sebelum", x = "Status", y = "Life Expectancy") +
  theme_minimal(base_size = 12)

p2 <- ggplot(data_plot, aes(x = Status, y = BMI, fill = IPM_Category)) +
  geom_boxplot(alpha = 0.8, position = position_dodge(width = 0.8), outlier.size = 0.8) +
  labs(title = "(b) BMI - Sebelum", x = "Status", y = "BMI") +
  theme_minimal(base_size = 12)

p3 <- ggplot(data_plot, aes(x = Status, y = Adult.Mortality, fill = IPM_Category)) +
  geom_boxplot(alpha = 0.8, position = position_dodge(width = 0.8), outlier.size = 0.8) +
  labs(title = "(c) Adult Mortality - Sebelum", x = "Status", y = "Adult Mortality") +
  theme_minimal(base_size = 12)

# Sesudah transformasi (YYY)
p4 <- ggplot(data_plot, aes(x = Status, y = Life_t, fill = IPM_Category)) +
  geom_boxplot(alpha = 0.8, position = position_dodge(width = 0.8), outlier.size = 0.8) +
  labs(title = "(d) Life Expectancy - Sesudah", x = "Status", y = "YYY$Life.expectancy") +
  theme_minimal(base_size = 12)

p5 <- ggplot(data_plot, aes(x = Status, y = BMI_t, fill = IPM_Category)) +
  geom_boxplot(alpha = 0.8, position = position_dodge(width = 0.8), outlier.size = 0.8) +
  labs(title = "(e) BMI - Sesudah", x = "Status", y = "YYY$BMI") +
  theme_minimal(base_size = 12)

p6 <- ggplot(data_plot, aes(x = Status, y = Mort_t, fill = IPM_Category)) +
  geom_boxplot(alpha = 0.8, position = position_dodge(width = 0.8), outlier.size = 0.8) +
  labs(title = "(f) Adult Mortality - Sesudah", x = "Status", y = "YYY$Adult.Mortality") +
  theme_minimal(base_size = 12)

grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3)

6.2 ASUMSI

6.2.1 UJI Normalitas

library(MVN)

# Sebelum transformasi / penanganan outlier
Y_before <- na.omit(
  data[, c("Life.expectancy", "BMI", "Adult.Mortality")]
)

result_before <- mvn(
  data = Y_before,
  mvn_test = "hz",
  univariate_test = "AD",
  tidy = FALSE
)

# Sesudah transformasi / penanganan outlier
Y_after <- na.omit(YYY)

result_after <- mvn(
  data = Y_after,
  mvn_test = "hz",
  univariate_test = "AD",
  tidy = FALSE
)

cat("<<<< Sebelum Transformasi <<<<\n")
## <<<< Sebelum Transformasi <<<<
result_before$multivariate_normality
##            Test Statistic p.value     Method
## 1 Henze-Zirkler  90.53056       0 asymptotic
result_before$univariate_normality
##               Test        Variable Statistic p.value
## 1 Anderson-Darling Life.expectancy  49.91453 3.7e-24
## 2 Anderson-Darling             BMI  78.50062 3.7e-24
## 3 Anderson-Darling Adult.Mortality  51.91747 3.7e-24
cat("\n<<<< Sesudah Transformasi <<<<\n")
## 
## <<<< Sesudah Transformasi <<<<
result_after$multivariate_normality
##            Test Statistic p.value     Method
## 1 Henze-Zirkler  90.53056       0 asymptotic
result_after$univariate_normality
##               Test        Variable Statistic p.value
## 1 Anderson-Darling Life.expectancy  49.91453 3.7e-24
## 2 Anderson-Darling             BMI  78.50062 3.7e-24
## 3 Anderson-Darling Adult.Mortality  51.91747 3.7e-24

6.2.2 Uji Independensi

library(psych)


# DATA ASLI
cor_before <- cor(
  data[, c("Life.expectancy", "BMI", "Adult.Mortality")],
  use = "complete.obs"
)

bart_before <- cortest.bartlett(
  cor_before,
  n = nrow(na.omit(data[, c("Life.expectancy", "BMI", "Adult.Mortality")]))
)

cat("<<<< Bartlett Test: Data Asli <<<<\n")
## <<<< Bartlett Test: Data Asli <<<<
print(bart_before)
## $chisq
## [1] 3064.846
## 
## $p.value
## [1] 0
## 
## $df
## [1] 3
# DATA TRANSFORMASI
cor_after <- cor(YYY, use = "complete.obs")

bart_after <- cortest.bartlett(
  cor_after,
  n = nrow(na.omit(YYY))
)

cat("\n<<<< Bartlett Test: Data Transformasi (YYY) <<<<\n")
## 
## <<<< Bartlett Test: Data Transformasi (YYY) <<<<
print(bart_after)
## $chisq
## [1] 3064.846
## 
## $p.value
## [1] 0
## 
## $df
## [1] 3

6.2.3 Uji Hommogenitas

library(biotools)


# DATA ASLI
Y_before <- na.omit(
  data[, c("Life.expectancy", "BMI", "Adult.Mortality")]
)

group_before <- data$Group[complete.cases(
  data[, c("Life.expectancy", "BMI", "Adult.Mortality")]
)]

cat("<<<< Box's M Test: Data Asli <<<<\n")
## <<<< Box's M Test: Data Asli <<<<
boxM(Y_before, group_before)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y_before
## Chi-Sq (approx.) = 2081.2, df = 42, p-value < 2.2e-16
# DATA TRANSFORMASI
Y_after <- na.omit(YYY)

group_after <- data$Group[complete.cases(YYY)]

cat("\n<<<< Box's M Test: Data Transformasi (YYY) <<<<\n")
## 
## <<<< Box's M Test: Data Transformasi (YYY) <<<<
boxM(Y_after, group_after)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y_after
## Chi-Sq (approx.) = 2081.2, df = 42, p-value < 2.2e-16

6.3 MODEL

6.3.1 Manova

# MANOVA TWO-WAY DATA ASLI

manova_before <- manova(
  cbind(Life.expectancy, BMI, Adult.Mortality) ~
    Status * IPM_Category,
  data = data
)

cat("<<<< MANOVA Two-Way: Data Asli <<<<\n")
## <<<< MANOVA Two-Way: Data Asli <<<<
summary(manova_before, test = "Pillai")
##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## Status                 1 0.44553   784.25      3   2928 < 2.2e-16 ***
## IPM_Category           3 0.62732   258.22      9   8790 < 2.2e-16 ***
## Status:IPM_Category    3 0.03993    13.18      9   8790 < 2.2e-16 ***
## Residuals           2930                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nANOVA univariat (Data Asli)\n")
## 
## ANOVA univariat (Data Asli)
summary.aov(manova_before)
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 2161.874 < 2.2e-16 ***
## IPM_Category           3 118282   39427 1385.602 < 2.2e-16 ***
## Status:IPM_Category    3   2592     864   30.363 < 2.2e-16 ***
## Residuals           2930  83373      28                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq F value    Pr(>F)    
## Status                 1 113688  113688 441.873 < 2.2e-16 ***
## IPM_Category           3 294160   98053 381.104 < 2.2e-16 ***
## Status:IPM_Category    3  10586    3529  13.715 6.985e-09 ***
## Residuals           2930 753851     257                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 419.5214 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 287.5568 < 2.2e-16 ***
## Status:IPM_Category    3   278144   92715   8.6813 9.849e-06 ***
## Residuals           2930 31291795   10680                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# MANOVA TWO-WAY DATA TRANSFORMASI (YYY)
data_after <- data
data_after$Life.expectancy <- YYY$Life.expectancy
data_after$BMI             <- YYY$BMI
data_after$Adult.Mortality <- YYY$Adult.Mortality

manova_after <- manova(
  cbind(Life.expectancy, BMI, Adult.Mortality) ~
    Status * IPM_Category,
  data = data_after
)

cat("\n<<<< MANOVA Two-Way: Data Transformasi (YYY) <<<<\n")
## 
## <<<< MANOVA Two-Way: Data Transformasi (YYY) <<<<
summary(manova_after, test = "Pillai")
##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## Status                 1 0.44553   784.25      3   2928 < 2.2e-16 ***
## IPM_Category           3 0.62732   258.22      9   8790 < 2.2e-16 ***
## Status:IPM_Category    3 0.03993    13.18      9   8790 < 2.2e-16 ***
## Residuals           2930                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\nANOVA univariat (Data Transformasi)\n")
## 
## ANOVA univariat (Data Transformasi)
summary.aov(manova_after)
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 2161.874 < 2.2e-16 ***
## IPM_Category           3 118282   39427 1385.602 < 2.2e-16 ***
## Status:IPM_Category    3   2592     864   30.363 < 2.2e-16 ***
## Residuals           2930  83373      28                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq F value    Pr(>F)    
## Status                 1 113688  113688 441.873 < 2.2e-16 ***
## IPM_Category           3 294160   98053 381.104 < 2.2e-16 ***
## Status:IPM_Category    3  10586    3529  13.715 6.985e-09 ***
## Residuals           2930 753851     257                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 419.5214 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 287.5568 < 2.2e-16 ***
## Status:IPM_Category    3   278144   92715   8.6813 9.849e-06 ***
## Residuals           2930 31291795   10680                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.3.2 Ancova

# ANCOVA DATA ASLI

ancova_before <- lm(
  Life.expectancy ~
    Status * IPM_Category +
    GDP + Schooling + HIV.AIDS,
  data = data
)

cat("<<<< ANCOVA: Data Asli <<<<\n")
## <<<< ANCOVA: Data Asli <<<<
anova(ancova_before)
## Analysis of Variance Table
## 
## Response: Life.expectancy
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANCOVA DATA TRANSFORMASI (YYY)

data_after <- data
data_after$Life.expectancy <- YYY$Life.expectancy

ancova_after <- lm(
  Life.expectancy ~
    Status * IPM_Category +
    GDP + Schooling + HIV.AIDS,
  data = data_after
)

cat("\n<<<< ANCOVA: Data Transformasi (YYY) <<<<\n")
## 
## <<<< ANCOVA: Data Transformasi (YYY) <<<<
anova(ancova_after)
## Analysis of Variance Table
## 
## Response: Life.expectancy
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.3.3 Mancova

# MANCOVA DATA ASLI
mancova_before <- manova(
  cbind(Life.expectancy, BMI, Adult.Mortality) ~
    Status * IPM_Category +
    GDP + Schooling + HIV.AIDS,
  data = data
)

cat("<<<< MANCOVA: Data Asli <<<<\n")
## <<<< MANCOVA: Data Asli <<<<
print(summary(mancova_before, test = "Pillai"))
##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## Status                 1 0.54928  1188.18      3   2925 < 2.2e-16 ***
## IPM_Category           3 0.71833   307.17      9   8781 < 2.2e-16 ***
## GDP                    1 0.00770     7.57      3   2925 4.843e-05 ***
## Schooling              1 0.06331    65.90      3   2925 < 2.2e-16 ***
## HIV.AIDS               1 0.36741   566.29      3   2925 < 2.2e-16 ***
## Status:IPM_Category    3 0.04490    14.83      9   8781 < 2.2e-16 ***
## Residuals           2927                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# MANCOVA DATA TRANSFORMASI
data_after <- data
data_after$Life.expectancy <- YYY$Life.expectancy
data_after$BMI             <- YYY$BMI
data_after$Adult.Mortality <- YYY$Adult.Mortality

mancova_after <- manova(
  cbind(Life.expectancy, BMI, Adult.Mortality) ~
    Status * IPM_Category +
    GDP + Schooling + HIV.AIDS,
  data = data_after
)

cat("\n<<<< MANCOVA: Data Transformasi (YYY) <<<<\n")
## 
## <<<< MANCOVA: Data Transformasi (YYY) <<<<
print(summary(mancova_after, test = "Pillai"))
##                       Df  Pillai approx F num Df den Df    Pr(>F)    
## Status                 1 0.54928  1188.18      3   2925 < 2.2e-16 ***
## IPM_Category           3 0.71833   307.17      9   8781 < 2.2e-16 ***
## GDP                    1 0.00770     7.57      3   2925 4.843e-05 ***
## Schooling              1 0.06331    65.90      3   2925 < 2.2e-16 ***
## HIV.AIDS               1 0.36741   566.29      3   2925 < 2.2e-16 ***
## Status:IPM_Category    3 0.04490    14.83      9   8781 < 2.2e-16 ***
## Residuals           2927                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("<<<< ANOVA Univariat dari MANCOVA Data Asli <<<<\n")
## <<<< ANOVA Univariat dari MANCOVA Data Asli <<<<
print(summary.aov(mancova_before))
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1 113688  113688 449.1656 < 2.2e-16 ***
## IPM_Category           3 294160   98053 387.3936 < 2.2e-16 ***
## GDP                    1    926     926   3.6576   0.05591 .  
## Schooling              1   7954    7954  31.4261 2.265e-08 ***
## HIV.AIDS               1   4144    4144  16.3739 5.334e-05 ***
## Status:IPM_Category    3  10560    3520  13.9068 5.293e-09 ***
## Residuals           2927 740853     253                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 520.9330 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 357.0684 < 2.2e-16 ***
## GDP                    1    27676   27676   3.2178   0.07294 .  
## Schooling              1     5007    5007   0.5822   0.44551    
## HIV.AIDS               1  6170617 6170617 717.4530 < 2.2e-16 ***
## Status:IPM_Category    3   192312   64104   7.4533 5.707e-05 ***
## Residuals           2927 25174326    8601                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n<<<< ANOVA Univariat dari MANCOVA Data Transformasi <<<<\n")
## 
## <<<< ANOVA Univariat dari MANCOVA Data Transformasi <<<<
print(summary.aov(mancova_after))
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1 113688  113688 449.1656 < 2.2e-16 ***
## IPM_Category           3 294160   98053 387.3936 < 2.2e-16 ***
## GDP                    1    926     926   3.6576   0.05591 .  
## Schooling              1   7954    7954  31.4261 2.265e-08 ***
## HIV.AIDS               1   4144    4144  16.3739 5.334e-05 ***
## Status:IPM_Category    3  10560    3520  13.9068 5.293e-09 ***
## Residuals           2927 740853     253                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 520.9330 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 357.0684 < 2.2e-16 ***
## GDP                    1    27676   27676   3.2178   0.07294 .  
## Schooling              1     5007    5007   0.5822   0.44551    
## HIV.AIDS               1  6170617 6170617 717.4530 < 2.2e-16 ***
## Status:IPM_Category    3   192312   64104   7.4533 5.707e-05 ***
## Residuals           2927 25174326    8601                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.3.4 Anova

# ANOVA UNIVARIAT SETELAH MANCOVA - DATA ASLI

cat("<<<< ANOVA Univariat MANCOVA: Data Asli <<<<\n")
## <<<< ANOVA Univariat MANCOVA: Data Asli <<<<
summary.aov(mancova_before)
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1 113688  113688 449.1656 < 2.2e-16 ***
## IPM_Category           3 294160   98053 387.3936 < 2.2e-16 ***
## GDP                    1    926     926   3.6576   0.05591 .  
## Schooling              1   7954    7954  31.4261 2.265e-08 ***
## HIV.AIDS               1   4144    4144  16.3739 5.334e-05 ***
## Status:IPM_Category    3  10560    3520  13.9068 5.293e-09 ***
## Residuals           2927 740853     253                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 520.9330 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 357.0684 < 2.2e-16 ***
## GDP                    1    27676   27676   3.2178   0.07294 .  
## Schooling              1     5007    5007   0.5822   0.44551    
## HIV.AIDS               1  6170617 6170617 717.4530 < 2.2e-16 ***
## Status:IPM_Category    3   192312   64104   7.4533 5.707e-05 ***
## Residuals           2927 25174326    8601                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA UNIVARIAT SETELAH MANCOVA - DATA TRANSFORMASI

cat("\n<<<< ANOVA Univariat MANCOVA: Data Transformasi (YYY) <<<<\n")
## 
## <<<< ANOVA Univariat MANCOVA: Data Transformasi (YYY) <<<<
summary.aov(mancova_after)
##  Response Life.expectancy :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  61516   61516 3379.719 < 2.2e-16 ***
## IPM_Category           3 118282   39427 2166.151 < 2.2e-16 ***
## GDP                    1    382     382   20.979 4.837e-06 ***
## Schooling              1   3014    3014  165.572 < 2.2e-16 ***
## HIV.AIDS               1  27395   27395 1505.063 < 2.2e-16 ***
## Status:IPM_Category    3   1899     633   34.779 < 2.2e-16 ***
## Residuals           2927  53276      18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response BMI :
##                       Df Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1 113688  113688 449.1656 < 2.2e-16 ***
## IPM_Category           3 294160   98053 387.3936 < 2.2e-16 ***
## GDP                    1    926     926   3.6576   0.05591 .  
## Schooling              1   7954    7954  31.4261 2.265e-08 ***
## HIV.AIDS               1   4144    4144  16.3739 5.334e-05 ***
## Status:IPM_Category    3  10560    3520  13.9068 5.293e-09 ***
## Residuals           2927 740853     253                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Adult.Mortality :
##                       Df   Sum Sq Mean Sq  F value    Pr(>F)    
## Status                 1  4480402 4480402 520.9330 < 2.2e-16 ***
## IPM_Category           3  9213143 3071048 357.0684 < 2.2e-16 ***
## GDP                    1    27676   27676   3.2178   0.07294 .  
## Schooling              1     5007    5007   0.5822   0.44551    
## HIV.AIDS               1  6170617 6170617 717.4530 < 2.2e-16 ***
## Status:IPM_Category    3   192312   64104   7.4533 5.707e-05 ***
## Residuals           2927 25174326    8601                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::

:::::::::::