library(readr)
## Warning: package 'readr' was built under R version 4.5.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(car)
## Warning: package 'car' was built under R version 4.5.3
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(MVN)
## Warning: package 'MVN' was built under R version 4.5.3
library(biotools)
## Warning: package 'biotools' was built under R version 4.5.3
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## ---
## biotools version 4.3
library(heplots)
## Warning: package 'heplots' was built under R version 4.5.3
## Loading required package: broom
## 
## Attaching package: 'heplots'
## The following object is masked from 'package:biotools':
## 
##     boxM
data2 <- read_csv("manufacturing_iot_quality_10000.csv")
## Rows: 10000 Columns: 20
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (6): plant, line_id, machine_id, shift, material_batch, defect_type
## dbl (14): record_id, operator_experience_years, ambient_temp_c, ambient_humi...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Lihat struktur data
str(data2)
## spc_tbl_ [10,000 × 20] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ record_id                : num [1:10000] 1 2 3 4 5 6 7 8 9 10 ...
##  $ plant                    : chr [1:10000] "Plant_A" "Plant_A" "Plant_A" "Plant_C" ...
##  $ line_id                  : chr [1:10000] "L3" "L1" "L4" "L3" ...
##  $ machine_id               : chr [1:10000] "M07" "M18" "M10" "M01" ...
##  $ shift                    : chr [1:10000] "evening" "night" "evening" "evening" ...
##  $ operator_experience_years: num [1:10000] 5.27 1.61 8.35 4.66 6.77 0 9.39 7.05 7.68 1.8 ...
##  $ material_batch           : chr [1:10000] "B090" "B122" "B071" "B069" ...
##  $ ambient_temp_c           : num [1:10000] 23.7 19 19.9 21.3 24.9 ...
##  $ ambient_humidity_pct     : num [1:10000] 67.9 34.8 50.4 34.8 32.7 90 48 55.1 41.9 44.7 ...
##  $ vibration_rms            : num [1:10000] 0.981 1.29 0.871 0.995 0.697 ...
##  $ pressure_bar             : num [1:10000] 5.08 4.8 5.64 5.28 5.05 ...
##  $ feed_rate_mm_s           : num [1:10000] 19.5 26.6 23.9 21.4 15 ...
##  $ laser_power_w            : num [1:10000] 304 377 329 323 341 ...
##  $ coolant_flow_l_min       : num [1:10000] 7.63 6.53 6.22 7.47 4.32 5.75 6.7 6.62 8.89 8.05 ...
##  $ thickness_mm             : num [1:10000] 1.22 1.28 1.22 1.2 1.22 ...
##  $ tensile_strength_mpa     : num [1:10000] 517 537 539 542 594 ...
##  $ surface_roughness_um     : num [1:10000] 2.68 2.81 2.14 2.2 1.34 ...
##  $ defect_flag              : num [1:10000] 0 0 0 0 0 0 0 1 0 0 ...
##  $ defect_type              : chr [1:10000] "none" "none" "none" "none" ...
##  $ rework_needed            : num [1:10000] 0 0 0 0 0 0 0 1 0 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   record_id = col_double(),
##   ..   plant = col_character(),
##   ..   line_id = col_character(),
##   ..   machine_id = col_character(),
##   ..   shift = col_character(),
##   ..   operator_experience_years = col_double(),
##   ..   material_batch = col_character(),
##   ..   ambient_temp_c = col_double(),
##   ..   ambient_humidity_pct = col_double(),
##   ..   vibration_rms = col_double(),
##   ..   pressure_bar = col_double(),
##   ..   feed_rate_mm_s = col_double(),
##   ..   laser_power_w = col_double(),
##   ..   coolant_flow_l_min = col_double(),
##   ..   thickness_mm = col_double(),
##   ..   tensile_strength_mpa = col_double(),
##   ..   surface_roughness_um = col_double(),
##   ..   defect_flag = col_double(),
##   ..   defect_type = col_character(),
##   ..   rework_needed = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>
head(data2)
## # A tibble: 6 × 20
##   record_id plant line_id machine_id shift operator_experience_…¹ material_batch
##       <dbl> <chr> <chr>   <chr>      <chr>                  <dbl> <chr>         
## 1         1 Plan… L3      M07        even…                   5.27 B090          
## 2         2 Plan… L1      M18        night                   1.61 B122          
## 3         3 Plan… L4      M10        even…                   8.35 B071          
## 4         4 Plan… L3      M01        even…                   4.66 B069          
## 5         5 Plan… L1      M08        even…                   6.77 B141          
## 6         6 Plan… L3      M13        day                     0    B049          
## # ℹ abbreviated name: ¹​operator_experience_years
## # ℹ 13 more variables: ambient_temp_c <dbl>, ambient_humidity_pct <dbl>,
## #   vibration_rms <dbl>, pressure_bar <dbl>, feed_rate_mm_s <dbl>,
## #   laser_power_w <dbl>, coolant_flow_l_min <dbl>, thickness_mm <dbl>,
## #   tensile_strength_mpa <dbl>, surface_roughness_um <dbl>, defect_flag <dbl>,
## #   defect_type <chr>, rework_needed <dbl>
colSums(is.na(data2))
##                 record_id                     plant                   line_id 
##                         0                         0                         0 
##                machine_id                     shift operator_experience_years 
##                         0                         0                         0 
##            material_batch            ambient_temp_c      ambient_humidity_pct 
##                         0                         0                         0 
##             vibration_rms              pressure_bar            feed_rate_mm_s 
##                         0                         0                         0 
##             laser_power_w        coolant_flow_l_min              thickness_mm 
##                         0                         0                         0 
##      tensile_strength_mpa      surface_roughness_um               defect_flag 
##                         0                         0                         0 
##               defect_type             rework_needed 
##                         0                         0
# Hapus missing value jika ada
data2_clean <- na.omit(data2)
data2_clean$shift <- as.factor(data2_clean$shift)

data2_clean$tensile_strength_mpa <- as.numeric(data2_clean$tensile_strength_mpa)
data2_clean$surface_roughness_um <- as.numeric(data2_clean$surface_roughness_um)
data2_clean$ambient_humidity_pct <- as.numeric(data2_clean$ambient_humidity_pct)
Y2 <- data2_clean[, c("tensile_strength_mpa",
                    "surface_roughness_um")]
# Uji Normalitas Univariat
# install.packages("nortest")
library(nortest)
ad.test(data2_clean$tensile_strength_mpa)
## 
##  Anderson-Darling normality test
## 
## data:  data2_clean$tensile_strength_mpa
## A = 0.27246, p-value = 0.6695
ad.test(data2_clean$surface_roughness_um)
## 
##  Anderson-Darling normality test
## 
## data:  data2_clean$surface_roughness_um
## A = 0.719, p-value = 0.0607
library(MVN)

Y2_clean <- data2_clean[, c("tensile_strength_mpa",
                            "surface_roughness_um")]

mvn_result <- mvn(
  data = Y2_clean,
  mvn_test = "mardia"
)

mvn_result$multivariate_normality
##              Test Statistic p.value     Method      MVN
## 1 Mardia Skewness     7.448   0.114 asymptotic ✓ Normal
## 2 Mardia Kurtosis     0.920   0.358 asymptotic ✓ Normal
# Homogenitas Matriks Kovarians
boxM(Y2, data2_clean$shift)
## 
##  Box's M-test for Homogeneity of Covariance Matrices 
## 
## data:  Y2 by data2_clean$shift 
## Chi-Sq (approx.) = 3.3034, df = 6, p-value = 0.7699
# Homogenitas Varians per DV
#--------------------------------------------------------
leveneTest(tensile_strength_mpa ~ shift, data = data2_clean)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    2  0.1823 0.8333
##       9997
leveneTest(surface_roughness_um ~ shift, data = data2_clean)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    2  1.8579  0.156
##       9997
plot(data2_clean$ambient_humidity_pct,
     data2_clean$tensile_strength_mpa)

plot(data2_clean$ambient_humidity_pct,
     data2_clean$surface_roughness_um)

model_mancova <- manova(
  cbind(tensile_strength_mpa,
        surface_roughness_um)
  ~ shift + ambient_humidity_pct,
  data = data2_clean
)
summary(model_mancova, test = "Pillai")
##                        Df   Pillai approx F num Df den Df    Pr(>F)    
## shift                   2 0.007486    18.78      4  19992 2.019e-15 ***
## ambient_humidity_pct    1 0.086195   471.39      2   9995 < 2.2e-16 ***
## Residuals            9996                                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Bisa juga:
summary(model_mancova, test = "Wilks")
##                        Df   Wilks approx F num Df den Df    Pr(>F)    
## shift                   2 0.99251    18.81      4  19990 1.894e-15 ***
## ambient_humidity_pct    1 0.91381   471.39      2   9995 < 2.2e-16 ***
## Residuals            9996                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(model_mancova)
##  Response tensile_strength_mpa :
##                        Df   Sum Sq Mean Sq  F value Pr(>F)    
## shift                   2      375     187   0.1744   0.84    
## ambient_humidity_pct    1   441559  441559 411.1148 <2e-16 ***
## Residuals            9996 10736240    1074                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response surface_roughness_um :
##                        Df  Sum Sq Mean Sq F value    Pr(>F)    
## shift                   2   12.49   6.245  37.523 < 2.2e-16 ***
## ambient_humidity_pct    1   89.45  89.455 537.521 < 2.2e-16 ***
## Residuals            9996 1663.54   0.166                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1