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
# load library
library(MVN)
library(biotools)
library(car)
# Dependen (Y)
y <- data[, c("Life.expectancy", "BMI", "Adult.Mortality")]
# Faktor
data$Status <- as.factor(data$Status)
# Kovariat (opsional, untuk nanti)
covars <- data[, c("GDP", "Schooling", "Income.composition.of.resources", "HIV.AIDS")]
result <- mvn(
data = y,
mvn_test = "hz",
univariate_test = "AD",
tidy = FALSE
)
result$multivariate_normality
## Test Statistic p.value Method
## 1 Henze-Zirkler 90.53056 0 asymptotic
result$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
#library(MVN)
#result <- MVN::mvn(data = y)
#result$multivariate_normality
model_final <- lm(
Life.expectancy ~ Status * (GDP + Schooling + Income.composition.of.resources + HIV.AIDS),
data = data
)
summary(model_final)
##
## Call:
## lm(formula = Life.expectancy ~ Status * (GDP + Schooling + Income.composition.of.resources +
## HIV.AIDS), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.4092 -2.5942 0.0371 2.6307 25.4536
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value
## (Intercept) 5.890e+01 2.689e+00 21.909
## StatusDeveloping -9.845e+00 2.716e+00 -3.625
## GDP 2.772e-05 1.074e-05 2.582
## Schooling -3.535e-02 1.624e-01 -0.218
## Income.composition.of.resources 2.426e+01 4.762e+00 5.095
## HIV.AIDS -7.184e-01 1.746e-02 -41.154
## StatusDeveloping:GDP 8.237e-05 1.652e-05 4.986
## StatusDeveloping:Schooling 1.278e+00 1.691e-01 7.562
## StatusDeveloping:Income.composition.of.resources -1.521e+01 4.817e+00 -3.158
## StatusDeveloping:HIV.AIDS NA NA NA
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## StatusDeveloping 0.000294 ***
## GDP 0.009877 **
## Schooling 0.827699
## Income.composition.of.resources 3.71e-07 ***
## HIV.AIDS < 2e-16 ***
## StatusDeveloping:GDP 6.53e-07 ***
## StatusDeveloping:Schooling 5.29e-14 ***
## StatusDeveloping:Income.composition.of.resources 0.001605 **
## StatusDeveloping:HIV.AIDS NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.638 on 2929 degrees of freedom
## Multiple R-squared: 0.7629, Adjusted R-squared: 0.7622
## F-statistic: 1178 on 8 and 2929 DF, p-value: < 2.2e-16
Y <- sqrt(y)
cor(Y)
## Life.expectancy BMI Adult.Mortality
## Life.expectancy 1.0000000 0.5073262 -0.6070671
## BMI 0.5073262 1.0000000 -0.3166680
## Adult.Mortality -0.6070671 -0.3166680 1.0000000
library(psych)
cor_mat <- cor(Y)
cortest.bartlett(cor_mat, n = nrow(data))
## $chisq
## [1] 2223.204
##
## $p.value
## [1] 0
##
## $df
## [1] 3
library(car)
leveneTest(Life.expectancy ~ Status, data = data)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 285.17 < 2.2e-16 ***
## 2936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(BMI ~ Status, data = data)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 240.13 < 2.2e-16 ***
## 2936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
leveneTest(Adult.Mortality ~ Status, data = data)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 261.38 < 2.2e-16 ***
## 2936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(biotools)
boxM(Y, data$Status)
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: Y
## Chi-Sq (approx.) = 652.83, df = 6, p-value < 2.2e-16
manova_model <- manova(
cbind(Life.expectancy, BMI, Adult.Mortality) ~ Status,
data = data
)
summary(manova_model, test = "Pillai")
## Df Pillai approx F num Df den Df Pr(>F)
## Status 1 0.23469 299.92 3 2934 < 2.2e-16 ***
## Residuals 2936
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(manova_model)
## Response Life.expectancy :
## Df Sum Sq Mean Sq F value Pr(>F)
## Status 1 61516 61516 884.28 < 2.2e-16 ***
## Residuals 2936 204247 70
## ---
## 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 315.31 < 2.2e-16 ***
## Residuals 2936 1058597 361
## ---
## 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 322.55 < 2.2e-16 ***
## Residuals 2936 40783082 13891
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ancova_model <- lm(
Life.expectancy ~ Status + GDP + Schooling,
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 1624.83 < 2.2e-16 ***
## GDP 1 17080 17080 451.14 < 2.2e-16 ***
## Schooling 1 76086 76086 2009.65 < 2.2e-16 ***
## Residuals 2934 111081 38
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mancova_model <- manova(
cbind(Life.expectancy, BMI, Adult.Mortality) ~
Status +
GDP + Schooling +
Income.composition.of.resources +
HIV.AIDS,
data = data
)
summary(mancova_model, test = "Pillai")
## Df Pillai approx F num Df den Df Pr(>F)
## Status 1 0.49221 946.69 3 2930 < 2.2e-16
## GDP 1 0.21215 263.00 3 2930 < 2.2e-16
## Schooling 1 0.55144 1200.67 3 2930 < 2.2e-16
## Income.composition.of.resources 1 0.12053 133.85 3 2930 < 2.2e-16
## HIV.AIDS 1 0.39003 624.50 3 2930 < 2.2e-16
## Residuals 2932
##
## Status ***
## GDP ***
## Schooling ***
## Income.composition.of.resources ***
## HIV.AIDS ***
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1