R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

library(psych)

## Warning: package 'psych' was built under R version 4.4.3

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.4.3

Input Data

berat_stall <- c(
  5.40, 5.14, 5.32, 5.28, 5.20, 5.26, 5.20, 5.24, 5.00, 5.00,
  5.56, 5.50, 5.20, 5.30, 5.20, 5.22, 5.40, 5.20, 5.16, 5.14,
  5.24, 5.10, 5.10, 5.20, 5.40, 5.30, 5.04, 5.20, 5.30, 5.20,
  5.06, 5.10, 5.16
)

tebal_stall <- c(
  0.92, 0.91, 0.91, 0.93, 0.90, 0.91, 0.90, 0.93, 0.92, 0.95,
  0.93, 0.90, 0.90, 0.90, 0.92, 0.90, 0.90, 0.93, 0.93, 0.93,
  0.93, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90, 0.90,
  0.87, 0.90, 0.90
)

df_stall <- data.frame(Berat = berat_stall, Tebal = tebal_stall)

berat_pipa <- c(
  3.80, 3.96, 3.72, 3.66, 3.96, 3.70, 3.80, 3.80, 3.72, 3.98,
  3.80, 3.80, 3.54, 3.74, 3.70, 3.61, 3.80, 3.80, 3.72, 3.72,
  3.80, 3.62, 3.66, 4.00, 3.60, 3.78, 3.80, 3.80, 3.80, 3.80,
  3.70, 3.68, 3.66, 3.90, 3.90
)

tebal_pipa <- c(
  1.15, 1.15, 1.15, 1.15, 1.10, 1.15, 1.10, 1.10, 1.15, 1.15,
  1.10, 1.10, 1.15, 1.15, 1.10, 1.15, 1.10, 1.10, 1.15, 1.10,
  1.10, 1.15, 1.15, 1.10, 0.90, 1.10, 1.10, 1.10, 1.10, 1.10,
  1.10, 1.10, 1.10, 1.10, 1.10
)

diameter_pipa <- c(
  24.5, 24.0, 24.0, 20.0, 24.0, 24.0, 24.0, 24.0, 21.0, 21.0,
  25.0, 24.0, 21.0, 21.0, 24.0, 21.0, 24.0, 24.0, 24.0, 24.0,
  24.0, 21.0, 21.0, 24.0, 24.0, 24.0, 22.0, 24.0, 24.0, 24.0,
  24.0, 21.0, 21.0, 25.0, 24.0
)

df_pipa <- data.frame(Berat = berat_pipa, Tebal = tebal_pipa, Diameter = diameter_pipa)

Statistika Deskriptif

desc_stats <- function(data) {
  data.frame(
    Mean     = colMeans(data),
    Varians  = apply(data, 2, var),
    Minimum  = apply(data, 2, min),
    Maksimum = apply(data, 2, max)
  )
}

cat("KARAKTERISTIK DATA STALL GALVANIS\n")

## KARAKTERISTIK DATA STALL GALVANIS

print(round(desc_stats(df_stall), 7))

##           Mean   Varians Minimum Maksimum
## Berat 5.221818 0.0171091    5.00     5.56
## Tebal 0.909697 0.0002530    0.87     0.95

cat("\nKARAKTERISTIK DATA PIPA GALVANIS\n")

## 
## KARAKTERISTIK DATA PIPA GALVANIS

print(round(desc_stats(df_pipa), 6))

##               Mean  Varians Minimum Maksimum
## Berat     3.766571 0.012158    3.54     4.00
## Tebal     1.112857 0.001962    0.90     1.15
## Diameter 23.128571 2.255042   20.00    25.00

Uji Asumsi Normalitas Multivariat (MARDIA TEST & QQ PLOT)

cat("\nUJI NORMALITAS MARDIA: STALL GALVANIS\n")

## 
## UJI NORMALITAS MARDIA: STALL GALVANIS

print(psych::mardia(df_stall, plot = FALSE))

## Call: psych::mardia(x = df_stall, plot = FALSE)
## 
## Mardia tests of multivariate skew and kurtosis
## Use describe(x) the to get univariate tests
## n.obs = 33   num.vars =  2 
## b1p =  1.25   skew =  6.85  with probability  <=  0.14
##  small sample skew =  7.94  with probability <=  0.094
## b2p =  9.04   kurtosis =  0.75  with probability <=  0.45

cat("\nUJI NORMALITAS MARDIA: PIPA GALVANIS\n")

## 
## UJI NORMALITAS MARDIA: PIPA GALVANIS

print(psych::mardia(df_pipa, plot = FALSE))

## Call: psych::mardia(x = df_pipa, plot = FALSE)
## 
## Mardia tests of multivariate skew and kurtosis
## Use describe(x) the to get univariate tests
## n.obs = 35   num.vars =  3 
## b1p =  14.07   skew =  82.05  with probability  <=  2e-13
##  small sample skew =  92.96  with probability <=  1.4e-15
## b2p =  27.28   kurtosis =  6.63  with probability <=  3.3e-11

# QQ Plot Mahalanobis Stall (df = 2)
d2_stall <- mahalanobis(df_stall, center = colMeans(df_stall), cov = cov(df_stall))
plot(sort(d2_stall), qchisq(ppoints(nrow(df_stall)), df = 2),
     xlab = "Squared Mahalanobis Distance", ylab = "X-Square Quantile",
     main = "QQ Plot Uji Normalitas Multivariat (Mahalanobis Distance)",
     pch = 19)
abline(0, 1, lwd = 2)
mtext("Derajat Kebebasan (df) = 2", side = 4, line = 0.5)

# QQ Plot Mahalanobis Pipa (df = 3)
d2_pipa <- mahalanobis(df_pipa, center = colMeans(df_pipa), cov = cov(df_pipa))
plot(qchisq(ppoints(nrow(df_pipa)), df = 3), sort(d2_pipa),
     xlab = "Kuantil Chi-Square", ylab = "Squared Mahalanobis Distance",
     main = "QQ Plot Uji Normalitas Multivariat (Mahalanobis Distance)",
     pch = 19)
abline(0, 1, lwd = 1)

# Uji Dependensi (BARTLETT SPHERICITY TEST)

cat("\nUJI BARTLETT: STALL GALVANIS\n")

## 
## UJI BARTLETT: STALL GALVANIS

print(cortest.bartlett(cor(df_stall), n = nrow(df_stall)))

## $chisq
## [1] 0.007133487
## 
## $p.value
## [1] 0.9326907
## 
## $df
## [1] 1

cat("\nUJI BARTLETT: PIPA GALVANIS\n")

## 
## UJI BARTLETT: PIPA GALVANIS

print(cortest.bartlett(cor(df_pipa), n = nrow(df_pipa)))

## $chisq
## [1] 14.6389
## 
## $p.value
## [1] 0.002152746
## 
## $df
## [1] 3

Peta Kendali Multivariat T^2 Hotelling

calc_t2_hotelling <- function(data, alpha = 0.001) {
  X <- as.matrix(data)
  n <- nrow(X)
  p <- ncol(X)
  X_bar <- colMeans(X)
  S_inv <- solve(cov(X))
  
  t2_values <- sapply(1:n, function(i) {
    diff <- X[i, ] - X_bar
    as.numeric(diff %*% S_inv %*% diff)
  })
  
  UCL <- ((p * (n - 1)) / (n - p)) * qf(1 - alpha, df1 = p, df2 = n - p)
  CL  <- mean(t2_values)
  
  return(list(T2 = t2_values, UCL = UCL, CL = CL))
} 

# Plot T^2 Stall
res_t2_stall <- calc_t2_hotelling(df_stall, alpha = 0.001)
plot(res_t2_stall$T2, type = "b", pch = 19, col = "blue",
     ylim = c(0, 17), xlab = "Sample", ylab = expression(T^2 ~ "Statistik"),
     main = expression("Peta Kendali " * T^2 * " Hotelling Stall"))
abline(h = 14, col = "red", lty = 2, lwd = 2)
abline(h = 1.9, col = "green3", lwd = 2)
text(x = 31, y = 15.5, labels = "UCL = 14", col = "red", font = 2)
text(x = 31, y = 2.7, labels = "CL = 1.9", col = "green3", font = 2)

# Plot T^2 Pipa 
res_t2_pipa <- calc_t2_hotelling(df_pipa, alpha = 0.001)
plot(res_t2_pipa$T2, type = "b", pch = 19, col = "navy",
     ylim = c(0, 27), xlab = "Sample", ylab = expression(T^2 ~ "Hotelling"),
     main = expression("Peta Kendali " * T^2 * " Hotelling Pipa\nVariabel: Berat, Tebal, Diameter. UCL pada alpha = 0.001"))
abline(h = round(res_t2_pipa$UCL, 2), col = "firebrick", lty = 2, lwd = 2)
abline(h = round(res_t2_pipa$CL, 2), col = "darkgreen", lty = 3, lwd = 2)
text(x = 32, y = 23, labels = paste0("UCL = ", round(res_t2_pipa$UCL, 2)), col = "firebrick")
text(x = 32, y = 3.8, labels = paste0("CL = ", round(res_t2_pipa$CL, 2)), col = "darkgreen")
points(25, res_t2_pipa$T2[25], col = "cyan3", pch = 19, cex = 1.5)
text(25, res_t2_pipa$T2[25] + 1.3, labels = "25", col = "red", font = 2)

# Stall Galvanis
w_lsl <- 5.0;  w_usl <- 5.5
t_lsl <- 0.9;  t_usl <- 1.1

mean_w_s <- mean(df_stall$Berat); sd_w_s <- sd(df_stall$Berat)
mean_t_s <- mean(df_stall$Tebal); sd_t_s <- sd(df_stall$Tebal)

cp_w_s  <- (w_usl - w_lsl) / (6 * sd_w_s)
cpk_w_s <- min((w_usl - mean_w_s) / (3 * sd_w_s), (mean_w_s - w_lsl) / (3 * sd_w_s))

cp_t_s  <- (t_usl - t_lsl) / (6 * sd_t_s)
cpk_t_s <- min((t_usl - mean_t_s) / (3 * sd_t_s), (mean_t_s - t_lsl) / (3 * sd_t_s))

mcp_stall  <- 0.5 * cp_w_s + 0.5 * cp_t_s
mcpk_stall <- 0.5 * cpk_w_s + 0.5 * cpk_t_s

cat("\nKAPABILITAS PROSES MULTIVARIAT STALL\n")

## 
## KAPABILITAS PROSES MULTIVARIAT STALL

print(data.frame(
  Variabel = c("Berat", "Tebal"),
  Wk       = c(0.5, 0.5),
  Cp       = round(c(cp_w_s, cp_t_s), 2),
  Cpk      = round(c(cpk_w_s, cpk_t_s), 2),
  MCp      = c(round(mcp_stall, 2), NA),
  MCpk     = c(round(mcpk_stall, 2), NA)
))

##   Variabel  Wk   Cp  Cpk  MCp MCpk
## 1    Berat 0.5 0.64 0.57 1.37 0.38
## 2    Tebal 0.5 2.10 0.20   NA   NA

# Pipa Galvanis 
pw_lsl <- 3.5;  pw_usl <- 4.0
pt_lsl <- 0.95; pt_usl <- 1.15
pd_lsl <- 21.0; pd_usl <- 25.0

mean_w_p <- mean(df_pipa$Berat);    sd_w_p <- sd(df_pipa$Berat)
mean_t_p <- mean(df_pipa$Tebal);    sd_t_p <- sd(df_pipa$Tebal)
mean_d_p <- mean(df_pipa$Diameter); sd_d_p <- sd(df_pipa$Diameter)

pp_w_p  <- (pw_usl - pw_lsl) / (6 * sd_w_p)
ppk_w_p <- min((pw_usl - mean_w_p) / (3 * sd_w_p), (mean_w_p - pw_lsl) / (3 * sd_w_p))

pp_t_p  <- (pt_usl - pt_lsl) / (6 * sd_t_p)
ppk_t_p <- min((pt_usl - mean_t_p) / (3 * sd_t_p), (mean_t_p - pt_lsl) / (3 * sd_t_p))

pp_d_p  <- (pd_usl - pd_lsl) / (6 * sd_d_p)
ppk_d_p <- min((pd_usl - mean_d_p) / (3 * sd_d_p), (mean_d_p - pd_lsl) / (3 * sd_d_p))

mpp_pipa  <- (1/3) * pp_w_p + (1/3) * pp_t_p + (1/3) * pp_d_p
mppk_pipa <- (1/3) * ppk_w_p + (1/3) * ppk_t_p + (1/3) * ppk_d_p

cat("\nKAPABILITAS PROSES MULTIVARIAT PIPA\n")

## 
## KAPABILITAS PROSES MULTIVARIAT PIPA

print(data.frame(
  Variabel = c("Berat", "Tebal", "Diameter"),
  Wk       = round(c(1/3, 1/3, 1/3), 3),
  Pp       = round(c(pp_w_p, pp_t_p, pp_d_p), 3),
  Ppk      = round(c(ppk_w_p, ppk_t_p, ppk_d_p), 3),
  MPp      = c(NA, round(mpp_pipa, 2), NA),
  MPpk     = c(NA, round(mppk_pipa, 2), NA)
))

##   Variabel    Wk    Pp   Ppk  MPp MPpk
## 1    Berat 0.333 0.756 0.706   NA   NA
## 2    Tebal 0.333 0.753 0.280 0.65 0.47
## 3 Diameter 0.333 0.444 0.415   NA   NA

VISUALISASI GGPLOT2 KAPABILITAS PROSES STALL

p1 <- ggplot(df_stall, aes(x = Berat)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.05, fill = "gray", color = "black") +
  stat_function(fun = dnorm, args = list(mean = mean_w_s, sd = sd_w_s),
                linewidth = 1, color = "black", linetype = "dashed") +
  geom_vline(xintercept = w_lsl, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = w_usl, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = w_lsl, y = 3.5, label = paste0("LSL: ", w_lsl), angle = 90, vjust = -0.5, color = "red", fontface = "bold") +
  annotate("text", x = w_usl, y = 3.5, label = paste0("USL: ", w_usl), angle = 90, vjust = 1.5, color = "red", fontface = "bold") +
  annotate("text", x = 5.55, y = 2.5, label = paste0("Cpk ≈ ", round(cpk_w_s, 2)), color = "black", hjust = 0) +
  labs(title = "Analisis Kapabilitas Proses - Berat (kg)",
       subtitle = paste0("Mean: ", round(mean_w_s, 4), " kg, StDev: ", round(sd_w_s, 4), " kg | Spesifikasi"),
       x = "Berat (kg)", y = "Density") +
  theme_minimal()

p2 <- ggplot(df_stall, aes(x = Tebal)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.01, fill = "gray", color = "black") +
  stat_function(fun = dnorm, args = list(mean = mean_t_s, sd = sd_t_s),
                linewidth = 1, color = "black", linetype = "dashed") +
  geom_vline(xintercept = t_lsl, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = t_usl, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = t_lsl, y = 30, label = paste0("LSL: ", t_lsl), angle = 90, vjust = -0.5, color = "red", fontface = "bold") +
  annotate("text", x = t_usl, y = 30, label = paste0("USL: ", t_usl), angle = 90, vjust = 1.5, color = "red", fontface = "bold") +
  annotate("text", x = 0.95, y = 30, label = paste0("Cpk ≈ ", round(cpk_t_s, 2)), color = "black", hjust = 0) +
  labs(title = "Analisis Kapabilitas Proses - Tebal (mm)",
       subtitle = paste0("(Mean: ", round(mean_t_s, 4), ", StDev: ", round(sd_t_s, 4), ")"),
       x = "Tebal (mm)", y = "Density") +
  theme_minimal()

grid.arrange(p1, p2, ncol = 2)

# Plot Kapabilitas - Berat Pipa
p_pipa_berat <- ggplot(df_pipa, aes(x = Berat)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.05, fill = "gray", color = "black") +
  stat_function(fun = dnorm, args = list(mean = mean_w_p, sd = sd_w_p),
                linewidth = 1, color = "black", linetype = "dashed") +
  geom_vline(xintercept = pw_lsl, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = pw_usl, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = pw_lsl, y = 1.0, label = paste0("LSL: ", pw_lsl), 
           angle = 90, vjust = -0.5, color = "red", fontface = "bold") +
  annotate("text", x = pw_usl, y = 1.0, label = paste0("USL: ", pw_usl), 
           angle = 90, vjust = 1.5, color = "red", fontface = "bold") +
  annotate("text", x = 5.5, y = 3.5, label = paste0("Cpk ≈ ", round(ppk_w_p, 2)), color = "black", hjust = 0) +
  xlim(3.5, 5.7) +
  labs(title = "Analisis Kapabilitas Proses - Berat (kg)",
       subtitle = paste0("Mean: ", round(mean_w_p, 4), " kg, StDev: ", round(sd_w_p, 4), " kg | Spesifikasi"),
       x = "Berat (kg)", y = "Density") +
  theme_minimal()

print(p_pipa_berat)

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

# Plot Kapabilitas - Tebal Pipa
p_pipa_tebal <- ggplot(df_pipa, aes(x = Tebal)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.02, fill = "gray", color = "black") +
  stat_function(fun = dnorm, args = list(mean = mean_t_p, sd = sd_t_p),
                linewidth = 1, color = "black", linetype = "dashed") +
  geom_vline(xintercept = pt_lsl, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = pt_usl, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = pt_lsl, y = 1.5, label = paste0("LSL: ", pt_lsl), 
           angle = 90, vjust = -0.5, color = "red", fontface = "bold") +
  annotate("text", x = pt_usl, y = 1.5, label = paste0("USL: ", pt_usl), 
           angle = 90, vjust = 1.5, color = "red", fontface = "bold") +
  annotate("text", x = 5.5, y = 3.5, label = paste0("Cpk ≈ ", round(ppk_t_p, 2)), color = "black", hjust = 0) +
  xlim(0.8, 5.7) +
  labs(title = "Analisis Kapabilitas Proses - Tebal (mm)",
       subtitle = paste0("Mean: ", round(mean_t_p, 4), " kg, StDev: ", round(sd_t_p, 4), " kg | Spesifikasi"),
       x = "Tebal (mm)", y = "Density") +
  theme_minimal()

print(p_pipa_tebal)

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

# Plot Kapabilitas - Diameter Pipa
p_pipa_diameter <- ggplot(df_pipa, aes(x = Diameter)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.5, fill = "gray", color = "black") +
  stat_function(fun = dnorm, args = list(mean = mean_d_p, sd = sd_d_p),
                linewidth = 1, color = "black", linetype = "dashed") +
  geom_vline(xintercept = pd_lsl, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = pd_usl, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = pd_lsl, y = 1.5, label = paste0("LSL: ", pd_lsl), 
           angle = 90, vjust = -0.5, color = "red", fontface = "bold") +
  annotate("text", x = pd_usl, y = 1.5, label = paste0("USL: ", pd_usl), 
           angle = 90, vjust = 1.5, color = "red", fontface = "bold") +
  annotate("text", x = 6.0, y = 3.5, label = "Cpk ≈ Ppk: 0.52", color = "black", hjust = 0) +
  xlim(5.5, 26) +
  labs(title = "Analisis Kapabilitas Proses - Diameter (mm)",
       subtitle = paste0("Mean: ", round(mean_d_p, 4), " kg, StDev: ", round(sd_d_p, 4), " kg | Spesifikasi"),
       x = "Diameter (mm)", y = "Density") +
  theme_minimal()

# Cetak grafik Diameter
print(p_pipa_diameter)

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).