UAS

Kontak : \(\downarrow\)
Email
Instagram https://www.instagram.com/jeremyheriyand/
RPubs https://rpubs.com/jeremyheriyandi23/
Nama Jeremy Saudi
NIM 20214920008
Prodi Statistika 2021

library(ggplot2)
library(dplyr)
library(broom)
library(ggpubr)
library(readxl)

Hypothesis Testing

Misalkan berat rata-rata Penguin Raja yang ditemukan di koloni Antartika tahun lalu adalah 15,4 kg. Dalam sampel 35 penguin pada waktu yang sama tahun ini di koloni yang sama, berat penguin rata-rata adalah 14,6 kg. Asumsikan standar deviasi populasi adalah 2,5 kg. Pada tingkat signifikansi 0,05, dapatkah kita menolak hipotesis nol bahwa berat penguin rata-rata tidak berbeda dari tahun lalu?

Hipotesis nol adalah bahwa μ = 15,4. Kita mulai dengan menghitung statistik pengujian.

xbar = 14.6            # sample mean 
mu0 = 15.4             # hypothesized value 
sigma = 2.5            # population standard deviation 
n = 35                 # sample size 
z = (xbar-mu0)/(sigma/sqrt(n)) 
z                      # test statisti
## [1] -1.893146

Maka hipotesis nol dari uji dua arah harus ditolak jika z≤−zα/2 atau z≥zα/2 .

Kami kemudian menghitung nilai kritis pada tingkat signifikansi 0,05.

alpha = .05 
z.half.alpha = qnorm(1-alpha/2) 
c(-z.half.alpha, z.half.alpha) 
## [1] -1.959964  1.959964

Statistik uji -1,8931 terletak di antara nilai kritis -1,9600 dan 1,9600. Oleh karena itu, pada tingkat signifikansi 0,05, kami tidak menolak hipotesis nol bahwa berat penguin rata-rata tidak berbeda dari tahun lalu.

pval = 2*pnorm(z)    # lower tail 
pval                   # two−tailed p−value 
## [1] 0.05833852

Karena ternyata lebih besar dari tingkat signifikansi 0,05, kami tidak menolak hipotesis nol bahwa μ = 15,4.

ANOVA

library(palmerpenguins)
library(tidyverse)

ping <- penguins %>%
  select(species, flipper_length_mm)

summary(ping)
##       species    flipper_length_mm
##  Adelie   :152   Min.   :172.0    
##  Chinstrap: 68   1st Qu.:190.0    
##  Gentoo   :124   Median :197.0    
##                  Mean   :200.9    
##                  3rd Qu.:213.0    
##                  Max.   :231.0    
##                  NA's   :2
boxplot(flipper_length_mm ~ species,
  data = ping
)

Ada satu outlier dalam kelompok Adelie, seperti yang didefinisikan oleh kriteria rentang interkuartil. Poin ini, bagaimanapun, tidak dilihat sebagai outlier yang signifikan sehingga kita dapat mengasumsikan bahwa asumsi tidak ada outlier signifikan terpenuhi.

# Boxplot
boxplot(flipper_length_mm ~ species,
  data = ping
)

Boxplot menunjukkan varians yang sama untuk spesies yang berbeda. Dalam boxplot, ini dapat dilihat dari fakta bahwa kotak dan kumis memiliki ukuran yang sebanding untuk semua spesies.

# Levene's test
library(car)

leveneTest(flipper_length_mm ~ species,
  data = ping
)

Nilai p lebih besar dari tingkat signifikansi 0,05, kita tidak menolak hipotesis nol, jadi kita tidak dapat menolak hipotesis bahwa varians sama antar spesies (nilai-p = 0,719).

Hasil ini juga sejalan dengan pendekatan visual, sehingga homogenitas varians terpenuhi baik secara visual maupun formal.

res_aov <- aov(flipper_length_mm ~ species,
  data = ping
)

summary(res_aov)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## species       2  52473   26237   594.8 <2e-16 ***
## Residuals   339  14953      44                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

Mengingat bahwa nilai-p lebih kecil dari 0,05, kami menolak hipotesis nol, jadi kami menolak hipotesis bahwa semua cara sama. Oleh karena itu, kita dapat menyimpulkan bahwa setidaknya satu spesies berbeda dari yang lain dalam hal panjang sirip (p-value < 2.2e-16).

(Demi ilustrasi, jika nilai-p lebih besar dari 0,05: kita tidak dapat menolak hipotesis nol bahwa semua cara adalah sama, jadi kita tidak dapat menolak hipotesis bahwa 3 spesies penguin yang dianggap sama dalam hal panjang sirip.)

MANOVA

inc = read.csv("prak6.csv",header=TRUE, sep=";")
summary (inc)
##     Panjang          Lebar            Tinggi           JK           
##  Min.   : 93.0   Min.   : 74.00   Min.   :35.00   Length:48         
##  1st Qu.:106.8   1st Qu.: 86.00   1st Qu.:40.00   Class :character  
##  Median :122.0   Median : 93.00   Median :44.50   Mode  :character  
##  Mean   :124.7   Mean   : 95.44   Mean   :46.38                     
##  3rd Qu.:136.5   3rd Qu.:102.00   3rd Qu.:51.00                     
##  Max.   :177.0   Max.   :132.00   Max.   :67.00
inc
# Summary statistics
inc %>%
  group_by(JK) %>%
  get_summary_stats(Panjang, Lebar, Tinggi, type = "mean_sd")
# Visualization
ggboxplot(
  inc, x = "JK", y = c("Panjang", "Lebar", "Tinggi"), 
  merge = TRUE, palette = "jco")+
  theme_minimal()

Check Sample Size Assumption

inc %>%
  group_by(JK) %>%
  summarise(N = n())

Karena tabel di atas menunjukkan 24 pengamatan per kelompok, asumsi ukuran sampel yang memadai terpenuhi.

Check Univariate Normality Assumption

library(rstatix)
inc %>%
  group_by(JK) %>%
  shapiro_test(Panjang, Lebar, Tinggi) %>%
  arrange(variable)

Jika data berdistribusi normal, nilai-p harus lebih besar dari 0,05. Karena data di atas p > 0.05, maka data berdistribusi normal.

Multivariate Normality

inc %>%
  select(Panjang, Lebar, Tinggi) %>%
  mshapiro_test()

The test is not significant (p > 0.05), so we can assume multivariate normality.

Identify Multicollinearity

Idealnya korelasi antara variabel hasil harus moderat, tidak terlalu tinggi. Korelasi di atas 0,9 merupakan indikasi multikolinearitas, yang bermasalah bagi MANOVA.

inc %>% cor_test(Panjang, Lebar, Tinggi)

Tidak ada multikolinearitas, seperti yang dinilai oleh korelasi Pearson ( < 0,0001). Dalam situasi di mana Anda memiliki multikolinearitas, Anda dapat mempertimbangkan untuk menghapus salah satu variabel hasil yang sangat berkorelasi.

Check Linearity Assumption

# Create a scatterplot matrix by group
library(GGally)
results <- inc %>%
  select(Panjang, Lebar, Tinggi, JK) %>%
  group_by(JK) %>%
  doo(~ggpairs(.) + theme_bw(), result = "plots")
results
results$plots
## [[1]]

## 
## [[2]]

Check the Homogeneity of Covariances Assumption

box_m(inc[, c("Panjang", "Lebar","Tinggi")], inc$JK)

The test is statistically significant (i.e., p < 0.001), so the data have violated the assumption of homogeneity of variance-covariance matrices.

Check the Homogneity of Variance Assumption

inc %>% 
  gather(key = "variable", value = "value", Panjang, Lebar, Tinggi) %>%
  group_by(variable) %>%
  levene_test(value ~ JK)

The Levene’s test is significant (p < 0.05), so there was no homogeneity of variances.

One-way MANOVA Test

#Compute MANOVA in an analysis of variance framework
Model1<- manova(cbind(Panjang, Lebar, Tinggi) ~ JK, inc)
summary(Model1)
##           Df  Pillai approx F num Df den Df    Pr(>F)    
## JK         1 0.61143   23.078      3     44 3.967e-09 ***
## Residuals 46                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ada perbedaan yang signifikan secara statistik antara JENIS KELAMIN pada variabel dependen gabungan , p < 0,0001. Dari output di atas, dapat dilihat bahwa ketiga variabel tersebut sangat berbeda secara signifikan antar Jenis Kelamin

Pairwise-comparison (One-way MANOVA)

#Compute MANOVA in a linear model framework
model2 <- lm(cbind(Panjang, Lebar, Tinggi) ~ JK, inc)
Manova(model2, test.statistic = "Pillai")
## 
## Type II MANOVA Tests: Pillai test statistic
##    Df test stat approx F num Df den Df    Pr(>F)    
## JK  1   0.61143   23.078      3     44 3.967e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Group the data by variable
grouped.data <- inc %>%
  gather(key = "variable", value = "value", Panjang, Lebar, Tinggi) %>%
  group_by(variable)

# Do welch one way anova test
grouped.data %>% welch_anova_test(value ~ JK)
# or do Kruskal-Wallis test
grouped.data %>% kruskal_test(value ~ JK)
# or use aov()
grouped.data %>% anova_test(value ~ JK)

From the results of anova_test(): There was a statistically significant difference , p < 0.0001 diantara Jenis kelamin pada kura-kura

Compute Multiple Pairwise Comparisons

pwc <- inc %>%
  gather(key = "variables", value = "value",Panjang, Lebar, Tinggi) %>%
  group_by(variables) %>%
  games_howell_test(value ~ JK) %>%
  select(-estimate, -conf.low, -conf.high) # Remove details
pwc

SEMUA NILAI PAIRWISE BERBEDA.