1 Persiapan Packages
library(lmtest)
library(sandwich)
library(tibble)
library(purrr)
library(ggplot2)
library(plotly)
library(tidyr)
2 Parameter dan Fungsi Simulasi
# Parameter eksplorasi
n_values <- c(1000, 2000, 3000, 4000, 5000)
sigma_values <- c(5, 10, 7, 15, 10)
r2_values <- c(0.5, 0.7, 0.9, 0.7, 0.9)
# Fungsi simulasi
simulasi_once <- function(n, sigma, r2) {
x1 <- rnorm(n, mean = 7, sd = 1)
x2 <- rnorm(n, mean = 3, sd = 1)
beta <- sqrt(r2 / (2 * (1 - r2)))
y <- beta * x1 + beta * x2 + rnorm(n, sd = sigma)
model1 <- lm(y ~ x1)
model2 <- lm(y ~ x2)
reset1 <- resettest(model1, power = 2:3, type = "regressor")
reset2 <- resettest(model2, power = 2:3, type = "regressor")
return(c(x1_sig = reset1$p.value < 0.05, x2_sig = reset2$p.value < 0.05))
}
power_test <- function(n, sigma, r2, reps = 500) {
hasil <- replicate(reps, simulasi_once(n, sigma, r2))
power_x1 <- mean(hasil[1, ])
power_x2 <- mean(hasil[2, ])
return(c(power_x1 = power_x1, power_x2 = power_x2))
}
3 Simulasi 125 kombinasi
kombinasi <- expand.grid(n = n_values, sigma = sigma_values, r2 = r2_values)
set.seed(123)
hasil_simulasi <- pmap_dfr(kombinasi, function(n, sigma, r2) {
power <- power_test(n, sigma, r2)
tibble(n = n, sigma = sigma, r2 = r2, power_x1 = power[1], power_x2 = power[2])
})
print(hasil_simulasi,n=125)
## # A tibble: 125 × 5
## n sigma r2 power_x1 power_x2
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1000 5 0.5 0.056 0.046
## 2 2000 5 0.5 0.048 0.034
## 3 3000 5 0.5 0.032 0.054
## 4 4000 5 0.5 0.052 0.048
## 5 5000 5 0.5 0.056 0.042
## 6 1000 10 0.5 0.052 0.028
## 7 2000 10 0.5 0.058 0.06
## 8 3000 10 0.5 0.068 0.068
## 9 4000 10 0.5 0.038 0.06
## 10 5000 10 0.5 0.042 0.062
## 11 1000 7 0.5 0.046 0.06
## 12 2000 7 0.5 0.06 0.042
## 13 3000 7 0.5 0.046 0.042
## 14 4000 7 0.5 0.052 0.072
## 15 5000 7 0.5 0.044 0.064
## 16 1000 15 0.5 0.066 0.04
## 17 2000 15 0.5 0.048 0.048
## 18 3000 15 0.5 0.066 0.06
## 19 4000 15 0.5 0.048 0.052
## 20 5000 15 0.5 0.044 0.054
## 21 1000 10 0.5 0.042 0.046
## 22 2000 10 0.5 0.064 0.038
## 23 3000 10 0.5 0.034 0.074
## 24 4000 10 0.5 0.044 0.07
## 25 5000 10 0.5 0.052 0.072
## 26 1000 5 0.7 0.066 0.058
## 27 2000 5 0.7 0.05 0.04
## 28 3000 5 0.7 0.03 0.036
## 29 4000 5 0.7 0.058 0.05
## 30 5000 5 0.7 0.056 0.048
## 31 1000 10 0.7 0.044 0.064
## 32 2000 10 0.7 0.064 0.052
## 33 3000 10 0.7 0.058 0.052
## 34 4000 10 0.7 0.044 0.04
## 35 5000 10 0.7 0.042 0.05
## 36 1000 7 0.7 0.062 0.05
## 37 2000 7 0.7 0.056 0.044
## 38 3000 7 0.7 0.05 0.046
## 39 4000 7 0.7 0.046 0.054
## 40 5000 7 0.7 0.06 0.028
## 41 1000 15 0.7 0.06 0.068
## 42 2000 15 0.7 0.062 0.054
## 43 3000 15 0.7 0.046 0.046
## 44 4000 15 0.7 0.052 0.05
## 45 5000 15 0.7 0.052 0.048
## 46 1000 10 0.7 0.054 0.07
## 47 2000 10 0.7 0.034 0.04
## 48 3000 10 0.7 0.046 0.05
## 49 4000 10 0.7 0.072 0.066
## 50 5000 10 0.7 0.052 0.064
## 51 1000 5 0.9 0.058 0.046
## 52 2000 5 0.9 0.044 0.042
## 53 3000 5 0.9 0.048 0.06
## 54 4000 5 0.9 0.042 0.042
## 55 5000 5 0.9 0.022 0.048
## 56 1000 10 0.9 0.036 0.052
## 57 2000 10 0.9 0.042 0.056
## 58 3000 10 0.9 0.06 0.052
## 59 4000 10 0.9 0.05 0.058
## 60 5000 10 0.9 0.038 0.068
## 61 1000 7 0.9 0.044 0.034
## 62 2000 7 0.9 0.052 0.056
## 63 3000 7 0.9 0.05 0.046
## 64 4000 7 0.9 0.046 0.034
## 65 5000 7 0.9 0.052 0.032
## 66 1000 15 0.9 0.048 0.042
## 67 2000 15 0.9 0.048 0.034
## 68 3000 15 0.9 0.07 0.048
## 69 4000 15 0.9 0.052 0.042
## 70 5000 15 0.9 0.054 0.044
## 71 1000 10 0.9 0.038 0.062
## 72 2000 10 0.9 0.04 0.042
## 73 3000 10 0.9 0.056 0.056
## 74 4000 10 0.9 0.038 0.034
## 75 5000 10 0.9 0.05 0.088
## 76 1000 5 0.7 0.058 0.046
## 77 2000 5 0.7 0.048 0.034
## 78 3000 5 0.7 0.05 0.048
## 79 4000 5 0.7 0.05 0.05
## 80 5000 5 0.7 0.032 0.058
## 81 1000 10 0.7 0.052 0.058
## 82 2000 10 0.7 0.06 0.05
## 83 3000 10 0.7 0.046 0.054
## 84 4000 10 0.7 0.06 0.04
## 85 5000 10 0.7 0.064 0.046
## 86 1000 7 0.7 0.062 0.048
## 87 2000 7 0.7 0.046 0.046
## 88 3000 7 0.7 0.054 0.034
## 89 4000 7 0.7 0.058 0.046
## 90 5000 7 0.7 0.044 0.042
## 91 1000 15 0.7 0.044 0.044
## 92 2000 15 0.7 0.042 0.06
## 93 3000 15 0.7 0.06 0.04
## 94 4000 15 0.7 0.044 0.042
## 95 5000 15 0.7 0.048 0.04
## 96 1000 10 0.7 0.052 0.078
## 97 2000 10 0.7 0.048 0.048
## 98 3000 10 0.7 0.048 0.048
## 99 4000 10 0.7 0.07 0.042
## 100 5000 10 0.7 0.046 0.064
## 101 1000 5 0.9 0.046 0.054
## 102 2000 5 0.9 0.058 0.048
## 103 3000 5 0.9 0.05 0.05
## 104 4000 5 0.9 0.046 0.052
## 105 5000 5 0.9 0.062 0.046
## 106 1000 10 0.9 0.032 0.056
## 107 2000 10 0.9 0.048 0.05
## 108 3000 10 0.9 0.046 0.05
## 109 4000 10 0.9 0.046 0.038
## 110 5000 10 0.9 0.054 0.05
## 111 1000 7 0.9 0.066 0.066
## 112 2000 7 0.9 0.05 0.034
## 113 3000 7 0.9 0.058 0.044
## 114 4000 7 0.9 0.044 0.058
## 115 5000 7 0.9 0.04 0.048
## 116 1000 15 0.9 0.06 0.052
## 117 2000 15 0.9 0.056 0.05
## 118 3000 15 0.9 0.048 0.068
## 119 4000 15 0.9 0.042 0.044
## 120 5000 15 0.9 0.06 0.048
## 121 1000 10 0.9 0.056 0.036
## 122 2000 10 0.9 0.044 0.04
## 123 3000 10 0.9 0.036 0.046
## 124 4000 10 0.9 0.064 0.06
## 125 5000 10 0.9 0.044 0.032
4 Visualisasi Heatmap
hasil_long <- hasil_simulasi %>%
pivot_longer(cols = c(power_x1, power_x2),
names_to = "variabel",
values_to = "power")
ggplot(hasil_long, aes(x = sigma, y = r2, fill = power)) +
geom_tile(color = "white") +
facet_grid(variabel ~ n) +
scale_fill_gradient(low = "white", high = "red") +
labs(title = "HeatMap Power Uji Ramsey RESET",
x = "Keragaman Galat (σ)",
y = "Keeratan Hubungan (R²)",
fill = "Power") +
theme_minimal()
5 Visualisasi 3D Plot
plot_ly(hasil_long,
x = ~sigma,
y = ~r2,
z = ~power,
color = ~variabel,
type = "scatter3d",
mode = "markers") %>%
layout(title = "Visualisasi 3D Power Uji Ramsey RESET",
scene = list(xaxis = list(title = "σ"),
yaxis = list(title = "R²"),
zaxis = list(title = "Power")))
6 Uji Ramsey RESET Kombinasi Manual dan Visualisasi Permukaan
(Opsional)
uji_kombinasi_manual_final <- function(n, sigma, r2, reps = 500, seed = 123) {
set.seed(seed)
beta <- sqrt(r2 / (2 * (1 - r2)))
simulasi_once <- function() {
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- beta * x1 + beta * x2 + rnorm(n, sd = sigma)
model1 <- lm(y ~ x1)
reset1 <- resettest(model1, power = 2:3, type = "regressor")
model2 <- lm(y ~ x2)
reset2 <- resettest(model2, power = 2:3, type = "regressor")
c(x1_sig = reset1$p.value < 0.05, x2_sig = reset2$p.value < 0.05)
}
hasil <- replicate(reps, simulasi_once())
power_x1 <- mean(hasil[1, ])
power_x2 <- mean(hasil[2, ])
cat("===== Hasil Power Uji Ramsey RESET (Manual) =====\n")
cat("Ukuran Sampel (n):", n, "\n")
cat("Standar Deviasi Galat (σ):", sigma, "\n")
cat("Keeratan Hubungan (R²):", r2, "\n")
cat("-------------------------------------------------\n")
cat("Power X1:", round(power_x1, 3), "\n")
cat("Power X2:", round(power_x2, 3), "\n")
# Ambil satu contoh data simulasi nyata
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- beta * x1 + beta * x2 + rnorm(n, sd = sigma)
model <- lm(y ~ x1 + x2)
# Visualisasi permukaan prediksi
x1_seq <- seq(-2, 2, length.out = 40)
x2_seq <- seq(-2, 2, length.out = 40)
grid <- expand.grid(x1 = x1_seq, x2 = x2_seq)
grid$y_hat <- predict(model, newdata = grid)
z_matrix <- matrix(grid$y_hat, nrow = length(x1_seq), ncol = length(x2_seq))
persp(x = x1_seq, y = x2_seq, z = z_matrix,
xlab = "x1", ylab = "x2", zlab = "ŷ",
main = paste0("Permukaan Prediksi (n=", n, ", σ=", sigma, ", R²=", r2, ")"),
theta = 30, phi = 20, col = "skyblue", shade = 0.5, border = NA)
}
uji_kombinasi_manual_final(n = 1000, sigma = 5, r2 = 0.5)
## ===== Hasil Power Uji Ramsey RESET (Manual) =====
## Ukuran Sampel (n): 1000
## Standar Deviasi Galat (σ): 5
## Keeratan Hubungan (R²): 0.5
## -------------------------------------------------
## Power X1: 0.056
## Power X2: 0.046
