# Set seed untuk reproduktibilitas
set.seed(123)
# Jumlah observasi
n <- 100
# Prediktor (X) dari distribusi normal
X <- rnorm(n, mean = 2, sd = 1)
# Parameter model
beta0 <- 0.5 # Intercept
beta1 <- 1.2 # Koefisien untuk X
# Mean (mu) dari model Poisson
mu <- exp(beta0 + beta1 * X)
# Menambahkan overdispersi dengan mengalikan varians
phi <- 2 # Parameter dispersi (overdispersi)
Y <- rpois(n, lambda = mu * phi)
# Gabungkan data ke dalam data frame
data <- data.frame(Y = Y, X = X)Quasi-Likelihood
Kasus 1: Data Count (Perhitungan Kejadian)
Berikut adalah simulasi untuk Kasus 1: Data Count menggunakan R. Kita akan membuat data count yang menunjukkan overdispersi dan memodelkannya menggunakan quasi-Poisson. Saya akan menjelaskan langkah demi langkahnya:
Langkah 1: Membuat Data Simulasi
Langkah 2: Mengecek Overdispersi
# Hitung mean dan varians dari Y
mean_Y <- mean(data$Y)
var_Y <- var(data$Y)
# Cek overdispersi
cat("Mean Y:", mean_Y, "\n")Mean Y: 73.6 cat("Varians Y:", var_Y, "\n")Varians Y: 9193.758 cat("Rasio Varians/Mean:", var_Y / mean_Y, "\n")Rasio Varians/Mean: 124.9152 Langkah 3: Fitting Model Quasi-Poisson
Kita akan menggunakan fungsi glm dengan keluarga quasipoisson untuk memodelkan data.
# Fitting model quasi-Poisson
model <- glm(Y ~ X, family = quasipoisson(link = "log"), data = data)
# Ringkasan model
summary(model)
Call:
glm(formula = Y ~ X, family = quasipoisson(link = "log"), data = data)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.19544 0.04229 28.27 <2e-16 ***
X 1.20185 0.01337 89.88 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 0.9223051)
Null deviance: 8359.074 on 99 degrees of freedom
Residual deviance: 89.771 on 98 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4Output yang perlu diperhatikan \(\beta\) : koeifisien untuk intercept dan prediktor, serta \(\phi\) : Dispersion parameter
berdasarkan output diatas didpat X berpengaruh terhadap log-mean dari Y, dan Dispersipn parameter 0.922 mendekati 1, artinya tidak terjadi overdispersi atau underdispersi.
Langkah 4: Prediksi dan Visualisasi
# Prediksi mean
data$predicted <- predict(model, type = "response")
# Plot data dan prediksi
plot(data$X, data$Y, pch = 16, col = "blue", xlab = "X", ylab = "Y", main = "Quasi-Poisson Model")
lines(data$X[order(data$X)], data$predicted[order(data$X)], col = "red", lwd = 2)
legend("topright", legend = c("Data", "Predicted"), col = c("blue", "red"), lty = c(NA, 1), pch = c(16, NA))
Langkah 5: Diagnosa Model
Periksa residu untuk memastikan kecocokan model.
# Hitung residu Pearson
residuals <- residuals(model, type = "pearson")
# Plot residu vs prediksi
plot(data$predicted, residuals, pch = 16, col = "blue", xlab = "Predicted", ylab = "Residuals", main = "Residuals vs Predicted")
abline(h = 0, col = "red", lty = 2)
Langkah 6: Estimasi Parameter Dispersi (\(\phi\) )
Parameter dispersi \(\phi\) dapat diestimasi secara manual.
# Estimasi parameter dispersi
phi_est <- sum(residuals^2) / (n - length(coef(model)))
cat("Estimated Dispersion Parameter (phi):", phi_est, "\n")Estimated Dispersion Parameter (phi): 0.9223051 Kasus 2: Data Proporsi
Berikut adalah simulasi untuk Kasus 2: Data Proporsi menggunakan R. Kita akan membuat data proporsi yang menunjukkan overdispersi dan memodelkannya menggunakan quasi-Binomial. Saya akan menjelaskan langkah demi langkahnya.
Langkah 1: Membuat Data Simulasi
# Set seed untuk reproduktibilitas
set.seed(123)
# Jumlah observasi
n <- 100
# Prediktor (X) dari distribusi normal
X <- rnorm(n, mean = 0, sd = 1)
# Parameter model
beta0 <- 0.5 # Intercept
beta1 <- 1.2 # Koefisien untuk X
# Probabilitas sukses (mu) dari model logit
mu <- plogis(beta0 + beta1 * X) # Fungsi logit invers
# Menambahkan overdispersi dengan mengalikan varians
phi <- 2 # Parameter dispersi (overdispersi)
# Menghasilkan data Binomial dengan overdispersi
size <- 10 # Jumlah percobaan untuk setiap observasi
Y <- rbinom(n, size = size, prob = mu * phi / (1 + mu * (phi - 1)))
# Ubah ke proporsi
Y_prop <- Y / size
# Gabungkan data ke dalam data frame
data <- data.frame(Y = Y_prop, X = X, Size = size)Langkah 2: Mengecek Overdispersi/Underdispersi
# Hitung mean dan varians dari Y
mean_Y <- mean(data$Y)
var_Y <- var(data$Y)
# Cek overdispersi
cat("Mean Y:", mean_Y, "\n")Mean Y: 0.745 cat("Varians Y:", var_Y, "\n")Varians Y: 0.05098485 cat("Rasio Varians/Mean:", var_Y / (mean_Y * (1 - mean_Y)), "\n")Rasio Varians/Mean: 0.2683766 Langkah 3: Fitting Model Quasi-Binomial
Kita akan menggunakan fungsi glm dengan keluarga quasibinomial untuk memodelkan data proporsi. Karena respons adalah proporsi, kita perlu menyertakan weights yang sesuai (dalam hal ini, jumlah percobaan size).
# Fitting model quasi-Binomial
model <- glm(Y ~ X, family = quasibinomial(link = "logit"), data = data, weights = Size)
# Ringkasan model
summary(model)
Call:
glm(formula = Y ~ X, family = quasibinomial(link = "logit"),
data = data, weights = Size)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.2535 0.0910 13.78 <2e-16 ***
X 1.3443 0.1182 11.38 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasibinomial family taken to be 1.083231)
Null deviance: 298.27 on 99 degrees of freedom
Residual deviance: 101.65 on 98 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 5Langkah 4: Prediksi dan Visualisasi
# Prediksi probabilitas
data$predicted <- predict(model, type = "response")
# Plot data dan prediksi
plot(data$X, data$Y, pch = 16, col = "blue", xlab = "X", ylab = "Y", main = "Quasi-Binomial Model for Proportion Data")
lines(data$X[order(data$X)], data$predicted[order(data$X)], col = "red", lwd = 2)
legend("topright", legend = c("Data", "Predicted"), col = c("blue", "red"), lty = c(NA, 1), pch = c(16, NA))
Langkah 5: Diagnosa Model
Periksa residu untuk memastikan kecocokan model.
# Hitung residu Pearson
residuals <- residuals(model, type = "pearson")
# Plot residu vs prediksi
plot(data$predicted, residuals, pch = 16, col = "blue", xlab = "Predicted", ylab = "Residuals", main = "Residuals vs Predicted")
abline(h = 0, col = "red", lty = 2)
Langkah 6: Estimasi Parameter Dispersi (\(\phi\))
Parameter dispersi \(\phi\) dapat diestimasi secara manual.
# Estimasi parameter dispersi
phi_est <- sum(residuals^2) / (n - length(coef(model)))
cat("Estimated Dispersion Parameter (phi):", phi_est, "\n")Estimated Dispersion Parameter (phi): 1.083231