Quasi Likelihood
Hutomo Atman Maulana, Yusma Yanti, Fahmi
STA2702 Inferensi Statistika Lanjut
A. Pendahuluan
Dalam Generalized Linear Model (GLM), pendekatan Maximum Likelihood Estimation (MLE) mengasumsikan bahwa data berasal dari keluarga eksponensial tertentu. Namun dalam praktik sering ditemukan fenomena , yaitu ketika varians empiris lebih besar daripada varians teoritis model.
Sebagai contoh: \[\begin{align} \text{Binomial:} \quad & \mathrm{Var}(Y_i) = n_i p_i (1-p_i) \\ \text{Poisson:} \quad & \mathrm{Var}(Y_i) = \mu_i \end{align}\]
Dalam banyak data nyata berlaku: \[\begin{equation} \mathrm{Var}(Y_i) > V(\mu_i) \end{equation}\]
Untuk mengatasi masalah ini, Wedderburn (1974) memperkenalkan pendekatan , yang hanya menspesifikasikan struktur mean dan varians tanpa memerlukan distribusi penuh.
B. Model Umum Quasi-Likelihood
Misalkan \(Y_1, \dots, Y_n\) adalah observasi independen.
Mean Model
\[\begin{equation} E(Y_i) = \mu_i \end{equation}\]
Variance Structure
\[\begin{equation} \mathrm{Var}(Y_i) = \phi V(\mu_i) \end{equation}\]
dengan:
Link Function
\[\begin{equation} g(\mu_i) = \eta_i = \mathbf{x}_i^T \mathbf{\beta} \end{equation}\]
Quasi-Likelihood
Quasi-likelihood \(Q(\mu_i; y_i)\) didefinisikan melalui turunan:
\[\begin{equation} \frac{\partial Q}{\partial \mu_i} = \frac{y_i - \mu_i}{V(\mu_i)} \end{equation}\]
Sehingga diperoleh:
\[\begin{equation} Q(\mu_i; y_i) = \int^{\mu_i} \frac{y_i - t}{V(t)} dt \end{equation}\]
Total quasi-likelihood:
\[\begin{equation} Q(\mathbf{\beta}) = \sum_{i=1}^n Q(\mu_i; y_i) \end{equation}\]
Quasi-Score Function
Turunan terhadap \(\mathbf{\beta}\) menghasilkan:
\[\begin{equation} U(\mathbf{\beta}) = \sum_{i=1}^n \frac{y_i - \mu_i}{V(\mu_i)} \frac{d\mu_i}{d\eta_i} \mathbf{x}_i \end{equation}\]
Dalam bentuk matriks:
\[\begin{equation} U(\mathbf{\beta}) = X^T D V^{-1} (\mathbf{y} - \mathbf{\mu}) \end{equation}\]
Fisher Scoring
Fisher Scoring:
\[\begin{equation} \mathbf{\beta}^{(t+1)} = \mathbf{\beta}^{(t)} + (X^T W X)^{-1} X^T W (\mathbf{y} - \mathbf{\mu}) \end{equation}\]
dengan bobot:
\[\begin{equation} W_i = \frac{1}{\phi} \frac{\left(\frac{d\mu_i}{d\eta_i}\right)^2} {V(\mu_i)} \end{equation}\]
Bentuk IRLS (Iteratively Reweighted Least Squares)
Definisikan :
\[\begin{equation} z_i = \eta_i + \frac{y_i - \mu_i} {\frac{d\mu_i}{d\eta_i}} \end{equation}\]
Update parameter:
\[\begin{equation} \mathbf{\beta}^{(t+1)} = (X^T W X)^{-1} X^T W \mathbf{z} \end{equation}\]
Estimasi Parameter Dispersi
Setelah konvergen, estimasi dispersi:
\[\begin{equation} \hat{\phi} = \frac{1}{n - p} \sum_{i=1}^n \frac{(y_i - \mu_i)^2} {V(\mu_i)} \end{equation}\]
Kovarians Estimator
\[\begin{equation} \mathrm{Var}(\hat{\mathbf{\beta}}) = \hat{\phi} (X^T W X)^{-1} \end{equation}\]
C. Algoritma Quasi-Likelihood
Misalkan tersedia data \(\{(y_i, \mathbf{x}_i)\}_{i=1}^n\) dengan model: \[ g(\mu_i) = \eta_i = \mathbf{x}_i^{T}\boldsymbol{\beta}, \qquad \mathrm{Var}(Y_i)=\phi V(\mu_i) \]
Tentukan:
Fungsi link\(g(\cdot)\)Fungsi varians\(V(\mu)\)
Inisialisasi Parameter:
Pilih nilai awal\(\boldsymbol{\beta}^{(0)}\).
Hitung Linear Predictor:\[ \eta_i^{(t)} = \mathbf{x}_i^{T}\boldsymbol{\beta}^{(t)} \]Hitung Mean Response:\[ \mu_i^{(t)} = g^{-1}(\eta_i^{(t)}) \]Hitung Turunan Mean:\[ \left( \frac{d\mu_i}{d\eta_i} \right)^{(t)} \]Hitung Bobot (Weight):\[ W_i^{(t)} = \frac{1}{\phi} \frac{\left(\frac{d\mu_i}{d\eta_i}\right)^2} {V(\mu_i)} \]
Bentuk matriks: \[
W^{(t)} = \mathrm{diag}(W_1^{(t)}, \dots, W_n^{(t)})
\]
Hitung Working Response:\[ z_i^{(t)} = \eta_i^{(t)} + \frac{y_i - \mu_i^{(t)}} {\left(\frac{d\mu_i}{d\eta_i}\right)^{(t)}} \]Update Parameter (IRLS):\[ \boldsymbol{\beta}^{(t+1)} = \left(X^T W^{(t)} X\right)^{-1} X^T W^{(t)} \mathbf{z}^{(t)} \]Cek Konvergensi:
Jika \[
\|\boldsymbol{\beta}^{(t+1)} - \boldsymbol{\beta}^{(t)}\| <
\varepsilon
\] maka berhenti;
jika tidak, kembali ke langkah 3.
Estimasi Parameter Dispersi Setelah konvergen:\[ \hat{\phi} = \frac{1}{n-p} \sum_{i=1}^{n} \frac{(y_i - \mu_i)^2} {V(\mu_i)} \]Kovarians Estimator:\[ \mathrm{Var}(\hat{\boldsymbol{\beta}}) = \hat{\phi} \left(X^T W X\right)^{-1} \]
Quasi-likelihood mempertahankan struktur mean model GLM namun memperlonggar asumsi distribusi penuh dengan menspesifikasikan fungsi varians dan parameter dispersi. Algoritma estimasi tetap menggunakan IRLS sebagaimana pada MLE, namun inferensi menjadi lebih robust terhadap overdispersion.
D. Simulasi
Data leaf blotch merupakan salah satu contoh klasik
dalam literatur Generalized Linear Models (GLM), yang diperkenalkan oleh
McCullagh dan Nelder (1989). Data ini memuat informasi mengenai proporsi
luas daun tanaman barley yang terinfeksi penyakit leaf blotch pada
berbagai kombinasi varietas dan lokasi tanam. Secara struktural, data
terdiri atas 90 observasi yang merepresentasikan hasil pengamatan pada
10 varietas barley yang diuji di 9 lokasi (site) yang
berbeda. Dengan demikian, struktur data bersifat faktorial dua arah, di
mana setiap kombinasi varietas dan lokasi menghasilkan satu pengukuran
proporsi infeksi.
Variabel respon dalam data ini adalah blotch, yang
menyatakan proporsi luas area daun yang terinfeksi. Nilai ini berada
dalam rentang \(0 < y < 1\).
Variabel penjelas terdiri atas dua faktor kategorik, yaitu
site (9 level) dan variety (10 level). Karena
respon berbentuk proporsi, model yang secara alami digunakan dalam
kerangka GLM adalah model binomial dengan fungsi link logit. Dalam
pendekatan ini diasumsikan bahwa
\[\begin{equation} Y_i \sim \text{Binomial}(n_i, p_i), \end{equation}\]
dengan model sistematik:
\[\begin{equation} \log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_{site} + \beta_{variety}. \end{equation}\]
Namun, dalam praktik sering ditemukan bahwa varians empiris data lebih besar daripada varians teoritis binomial, yaitu
\[\begin{equation} \mathrm{Var}(Y_i) > n_i p_i(1-p_i). \end{equation}\]
Fenomena ini dikenal sebagai . Untuk mengakomodasi kondisi tersebut, digunakan pendekatan quasi-binomial yang memperkenalkan parameter dispersi \(\phi\), sehingga struktur varians menjadi
\[\begin{equation} \mathrm{Var}(Y_i) = \phi n_i p_i(1-p_i). \end{equation}\]
Model quasi-binomial mempertahankan struktur mean dan link function yang sama seperti model binomial, tetapi memungkinkan varians lebih fleksibel melalui estimasi parameter dispersi dari data. Dengan demikian, analisis dapat menghasilkan standar error yang lebih realistis dan inferensi yang lebih robust terhadap pelanggaran asumsi distribusi binomial murni. Secara keseluruhan, data leaf blotch memberikan contoh yang sangat baik untuk mengilustrasikan perbedaan antara pendekatan binomial klasik dan quasi-likelihood, terutama dalam konteks overdispersion dan penyesuaian inferensi statistik.
# Load libraries
library(faraway)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(gnm)##
## Attaching package: 'gnm'## The following object is masked from 'package:faraway':
##
## wheat# Load the data leafblotch
data(leafblotch)
glimpse(leafblotch)## Rows: 90
## Columns: 3
## $ blotch <dbl> 0.0005, 0.0000, 0.0000, 0.0010, 0.0025, 0.0005, 0.0050, 0.0130…
## $ site <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,…
## $ variety <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, …#Explore data
summary(leafblotch)## blotch site variety
## Min. :0.0000 1 :10 1 : 9
## 1st Qu.:0.0125 2 :10 2 : 9
## Median :0.0500 3 :10 3 : 9
## Mean :0.2072 4 :10 4 : 9
## 3rd Qu.:0.3550 5 :10 5 : 9
## Max. :0.9500 6 :10 6 : 9
## (Other):30 (Other):36#Put the data in table
DataTable=matrix(leafblotch$blotch,ncol=10,
dimnames=list(paste("site",1:9,sep=""),paste("vrty",1:10,sep="")))
DataTable## vrty1 vrty2 vrty3 vrty4 vrty5 vrty6 vrty7 vrty8 vrty9 vrty10
## site1 0.0005 0.0150 0.0100 0.2000 0.250 0.0800 0.050 0.05 0.050 0.250
## site2 0.0000 0.0000 0.1270 0.3750 0.550 0.1658 0.050 0.50 0.050 0.425
## site3 0.0000 0.0005 0.0125 0.2625 0.050 0.2950 0.100 0.10 0.250 0.500
## site4 0.0010 0.0005 0.0125 0.0250 0.400 0.2000 0.050 0.50 0.750 0.375
## site5 0.0025 0.0030 0.0250 0.0050 0.055 0.4350 0.500 0.25 0.500 0.950
## site6 0.0005 0.0075 0.1660 0.0001 0.010 0.0100 0.750 0.50 0.750 0.625
## site7 0.0050 0.0030 0.0250 0.0300 0.060 0.0500 0.050 0.75 0.750 0.950
## site8 0.0130 0.0300 0.0250 0.0250 0.011 0.0500 0.001 0.05 0.750 0.950
## site9 0.0150 0.0750 0.0000 0.0001 0.025 0.0500 0.050 0.10 0.175 0.950# Define variable
blotch <- leafblotch$blotch
site <- leafblotch$site
variety <- leafblotch$varietyPemodelan dengan Binomial GLM
# Define n (untuk Binomial)
n<-1000
# Fit model GLM dengan family binomial dengan weights=n
model_bin <- glm(blotch ~ site + variety, data = leafblotch,
family = binomial, weights = rep(n, nrow(leafblotch)))## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(model_bin)##
## Call:
## glm(formula = blotch ~ site + variety, family = binomial, data = leafblotch,
## weights = rep(n, nrow(leafblotch)))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.05460 0.15091 -53.372 <2e-16 ***
## site2 1.63906 0.15318 10.700 <2e-16 ***
## site3 3.32648 0.14320 23.230 <2e-16 ***
## site4 3.58219 0.14269 25.105 <2e-16 ***
## site5 3.58384 0.14268 25.117 <2e-16 ***
## site6 3.89320 0.14224 27.371 <2e-16 ***
## site7 4.72994 0.14166 33.389 <2e-16 ***
## site8 5.52259 0.14165 38.989 <2e-16 ***
## site9 6.79450 0.14229 47.751 <2e-16 ***
## variety2 0.15008 0.07681 1.954 0.0507 .
## variety3 0.68946 0.07136 9.662 <2e-16 ***
## variety4 1.04814 0.06893 15.207 <2e-16 ***
## variety5 1.61468 0.06640 24.316 <2e-16 ***
## variety6 2.37110 0.06464 36.684 <2e-16 ***
## variety7 2.57117 0.06437 39.946 <2e-16 ***
## variety8 3.34193 0.06384 52.347 <2e-16 ***
## variety9 3.49994 0.06382 54.839 <2e-16 ***
## variety10 4.25287 0.06413 66.319 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 40802.9 on 89 degrees of freedom
## Residual deviance: 6126.4 on 72 degrees of freedom
## AIC: 6637.7
##
## Number of Fisher Scoring iterations: 5# Cek residuals untuk model Binomial
res_pearson <- residuals(model_bin, type="pearson")
linpred <- predict(model_bin, type="link")
# Plot 4 Panel
par(mfrow=c(2,2), mar=c(4,4,3,1))
### (1) Residual vs Linear Predictor
plot(linpred, res_pearson,
main="Pearson Residuals",
xlab="Linear Predictor",
ylab="Residual",
col="blue", pch=19)
abline(h=0, col="gray", lwd=2)
### (2) Histogram + Density
hist(res_pearson,
probability=TRUE,
main="",
xlab="Residual",
col="lightgray",
border="white")
lines(density(res_pearson),
col="blue", lwd=2)
### (3) QQ-Plot
qqnorm(res_pearson,
main="")
qqline(res_pearson,
col="gray", lwd=2)
### (4) Boxplot
boxplot(res_pearson,
horizontal=FALSE,
main="",
col="lightgray")
abline(h=0, col="gray", lwd=2)
# Reset layout
par(mfrow=c(1,1))plot(model_bin,which=1)
# Cek overdispersion dengan pearson residual deviance
pearson_residuals <- residuals(model_bin, type = "pearson")
deviance <- sum(pearson_residuals^2)
degrees_of_freedom <- df.residual(model_bin)
overdispersion <- deviance / degrees_of_freedom
overdispersion## [1] 88.78094Pemodelan dengan Quasi-Binomial GLM
# Jika overdispersion > 1, maka menggunakan model quasibinomial standar
model_quasi <- glm(blotch ~ site+variety,
family = quasi(link = "logit", variance = "mu(1-mu)"),
data = leafblotch)
summary(model_quasi)##
## Call:
## glm(formula = blotch ~ site + variety, family = quasi(link = "logit",
## variance = "mu(1-mu)"), data = leafblotch)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.0546 1.4220 -5.664 2.84e-07 ***
## site2 1.6391 1.4433 1.136 0.259881
## site3 3.3265 1.3492 2.465 0.016069 *
## site4 3.5822 1.3445 2.664 0.009512 **
## site5 3.5838 1.3444 2.666 0.009479 **
## site6 3.8932 1.3402 2.905 0.004877 **
## site7 4.7299 1.3348 3.544 0.000698 ***
## site8 5.5226 1.3346 4.138 9.39e-05 ***
## site9 6.7945 1.3407 5.068 3.00e-06 ***
## variety2 0.1501 0.7237 0.207 0.836293
## variety3 0.6895 0.6724 1.025 0.308599
## variety4 1.0481 0.6494 1.614 0.110919
## variety5 1.6147 0.6257 2.581 0.011897 *
## variety6 2.3711 0.6090 3.893 0.000219 ***
## variety7 2.5712 0.6065 4.240 6.55e-05 ***
## variety8 3.3419 0.6015 5.556 4.39e-07 ***
## variety9 3.4999 0.6014 5.820 1.51e-07 ***
## variety10 4.2529 0.6042 7.038 9.39e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasi family taken to be 0.08878102)
##
## Null deviance: 40.8029 on 89 degrees of freedom
## Residual deviance: 6.1264 on 72 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5# Cek residuals untuk model Quasi
res_pearson_quasi <- residuals(model_quasi, type="pearson")
linpred_quasi <- predict(model_quasi, type="link")
# Plot 4 Panel
par(mfrow=c(2,2), mar=c(4,4,3,
1))
### (1) Residual vs Linear Predictor
plot(linpred_quasi, res_pearson_quasi,
main="Pearson Residuals (Quasi)",
xlab="Linear Predictor",
ylab="Residual",
col="blue", pch=19)
abline(h=0, col="gray", lwd=2)
### (2) Histogram + Density
hist(res_pearson_quasi,
probability=TRUE,
main="",
xlab="Residual",
col="lightgray",
border="white")
lines(density(res_pearson_quasi),
col="blue", lwd=2)
### (3) QQ-Plot
qqnorm(res_pearson_quasi,
main="")
qqline(res_pearson_quasi,
col="gray", lwd=2)
### (4) Boxplot
boxplot(res_pearson_quasi,
horizontal=FALSE,
main="",
xlab="Residual",
col="lightgray")
abline(h=0, col="gray", lwd=2)
#Reset layout
par(mfrow=c(1,1))plot(model_quasi,which=1)
# Menghitung Quasi AIC
# Dispersion
phi <- summary(model_quasi)$dispersion
# Deviance
dev <- model_quasi$deviance
# Jumlah parameter
k <- length(coef(model_quasi))
# QAIC
QAIC <- dev/phi + 2*k
QAIC## [1] 105.0062#membandingkan QAIC dengan AIC dari model
AIC_bin <- AIC(model_bin)
AIC_quasi <- QAIC
cat("AIC (Binomial):", AIC_bin, "\n")## AIC (Binomial): 6637.671cat("QAIC (Quasi):", AIC_quasi, "\n")## QAIC (Quasi): 105.0062# Perbandingan Chi-Square Test
chi_bin=sum(resid(model_bin,type="pearson")^2)
chi_quasi=sum(resid(model_quasi,type="pearson")^2)
#buat tabel perbandingan
comparison_table <- data.frame(
Model = c("Binomial", "Quasi"),
Chi_Square = c(chi_bin, chi_quasi)
)
print(comparison_table)## Model Chi_Square
## 1 Binomial 6392.227933
## 2 Quasi 6.392228Custom Variance Function
#Custom Variance Function
CustomVar <- list(
varfun = function(mu) (mu * (1 - mu))^2,
validmu = function(mu) all(mu > 0) && all(mu < 1),
dev.resids = function(y, mu, wt) wt * ((y - mu)^2) / (mu * (1 - mu))^2,
initialize = expression({
n <- rep.int(1, nobs)
mustart <- pmax(0.001, pmin(0.999, y))
}),
name = "(mu(1-mu))^2")
# Fit model dengan custom variance function
model_custom <- glm(blotch ~ site + variety,
family = quasi(link = "logit", variance = CustomVar),
data = leafblotch)
summary(model_custom)##
## Call:
## glm(formula = blotch ~ site + variety, family = quasi(link = "logit",
## variance = CustomVar), data = leafblotch)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.92253 0.44463 -17.818 < 2e-16 ***
## site2 1.38308 0.44463 3.111 0.00268 **
## site3 3.86013 0.44463 8.682 8.18e-13 ***
## site4 3.55697 0.44463 8.000 1.53e-11 ***
## site5 4.10841 0.44463 9.240 7.48e-14 ***
## site6 4.30541 0.44463 9.683 1.13e-14 ***
## site7 4.91811 0.44463 11.061 < 2e-16 ***
## site8 5.69492 0.44463 12.808 < 2e-16 ***
## site9 7.06762 0.44463 15.896 < 2e-16 ***
## variety2 -0.46728 0.46868 -0.997 0.32210
## variety3 0.07877 0.46868 0.168 0.86699
## variety4 0.95419 0.46868 2.036 0.04544 *
## variety5 1.35276 0.46868 2.886 0.00514 **
## variety6 1.32859 0.46868 2.835 0.00595 **
## variety7 2.34066 0.46868 4.994 3.99e-06 ***
## variety8 3.26268 0.46868 6.961 1.30e-09 ***
## variety9 3.13556 0.46868 6.690 4.10e-09 ***
## variety10 3.88736 0.46868 8.294 4.33e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasi family taken to be 0.9884757)
##
## Null deviance: 252.14 on 89 degrees of freedom
## Residual deviance: 71.17 on 72 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 20#Menghitung dispersi untuk model custom
phi_custom <- summary(model_custom)$dispersion
cat("Estimated Dispersion (Custom):", phi_custom, "\n")## Estimated Dispersion (Custom): 0.9884757#Plot residuals untuk model custom
plot(model_custom,which=1) #Deviance resids vs linear predictors
# Menghitung QAIC untuk model custom
phi_c <- summary(model_custom)$dispersion
# Deviance
dev_c <- model_custom$deviance
# Jumlah parameter
k_c <- length(coef(model_custom))
# QAIC
QAIC_custom <- dev_c/phi_c + 2*k_c
QAIC_custom## [1] 108# Perbandingan model custom dengan model sebelumnya
comparison_table_c <- data.frame(
Model = c("Binomial", "Quasi_Binomial", "Quasi_Custom"),
Dispersion = c(1, phi, phi_c),
Variance_Function = c("mu(1-mu)", "mu(1-mu)", "(mu(1-mu))^2"),
AIC_QAIC = c(AIC_bin, QAIC, QAIC_custom),
Chi_Square = c(chi_bin, chi_quasi, sum(resid(model_custom,type="pearson")^2))
)
print(comparison_table_c)## Model Dispersion Variance_Function AIC_QAIC Chi_Square
## 1 Binomial 1.00000000 mu(1-mu) 6637.6713 6392.227933
## 2 Quasi_Binomial 0.08878102 mu(1-mu) 105.0062 6.392228
## 3 Quasi_Custom 0.98847570 (mu(1-mu))^2 108.0000 71.170250E. Interpretasi Hasil
Ringkasan Perbandingan Model
Tiga model yang dibandingkan adalah:
Ringkasan hasil estimasi:
Karena model quasi tidak memiliki likelihood penuh, maka pembanding yang tepat adalah QAIC. Model terbaik dipilih berdasarkan nilai QAIC terkecil.
Interpretasi Berdasarkan Struktur Model
1. Model Binomial (MLE)
Model Binomial mengasumsikan dispersion = 1 dengan struktur varians:
\[ Var(Y_i) = n_i p_i (1-p_i). \]
Nilai AIC yang sangat besar (6637.6713) serta nilai Chi-Square yang tinggi (6392.227933) menunjukkan bahwa asumsi varians binomial tidak sesuai dengan data.
Konsekuensi empiris:
2. Model Quasi-Binomial
Model ini mempertahankan struktur mean dan variance function \(\mu(1-\mu)\), namun mengestimasi parameter dispersi dari data.
\[ Var(Y_i) = \phi n_i p_i (1-p_i) \]
Nilai dispersion = 0.08878102 menunjukkan adanya penyesuaian varians terhadap data. QAIC yang paling kecil (105.0062) dan Chi-Square yang sangat kecil (6.392228) menunjukkan model ini paling sesuai.
Konsekuensi empiris:
3. Model Quasi-Custom Variance
Model quasi-custom menggunakan struktur fungsi varians berbentuk kuadrat,
\[ V(\mu) = (\mu(1-\mu))^2, \]
sehingga model lengkapnya menjadi
\[ Var(Y_i) = \phi \, (\mu_i(1-\mu_i))^2. \]
Struktur ini berbeda dari model quasi-binomial standar yang menggunakan \(V(\mu)=\mu(1-\mu)\), karena pada model ini varians meningkat secara lebih tajam terhadap perubahan mean. Secara empiris, bentuk kuadrat tersebut lebih mampu mengakomodasi pola variabilitas data leaf blotch yang tidak mengikuti struktur binomial klasik.
Hasil estimasi menunjukkan bahwa nilai dispersion mendekati satu (\(\hat{\phi} \approx 0.989\)), serta residual deviance sebesar 71.175 sangat dekat dengan derajat bebas model (72). Kondisi ini mengindikasikan bahwa struktur varians yang diasumsikan telah menangkap variasi data secara memadai.
Meskipun nilai QAIC model quasi-custom (108.0000) sedikit lebih besar dibanding quasi-binomial (105.0062), pemilihan model dalam kerangka quasi-likelihood tidak semata-mata ditentukan oleh kriteria informasi. Pertimbangan utama adalah kesesuaian struktur varians terhadap karakteristik empiris data.
Dalam konteks ini, model quasi-custom lebih konsisten secara teoretis dan metodologis, sebagaimana dijelaskan oleh Millar (2011), karena memberikan representasi varians yang lebih realistis dan menghasilkan inferensi yang lebih stabil terhadap pola heterogenitas dalam data.
F. Kesimpulan
Berdasarkan hasil estimasi tiga model, yaitu: Binomial GLM, Quasi-Binomial dengan \(V(\mu)=\mu(1-\mu)\), dan Quasi-Custom dengan \(V(\mu)=(\mu(1-\mu))^2\) terdapat perbedaan penting baik secara statistik maupun konseptual.
Model binomial menghasilkan nilai AIC yang sangat besar dan nilai Pearson chi-square yang tinggi, yang menunjukkan bahwa asumsi varians binomial klasik tidak sesuai dengan data. Hal ini mengindikasikan adanya pada struktur varians, sehingga inferensi dari model binomial cenderung terlalu optimistik dan standard error menjadi tidak realistis.
Model quasi-binomial memberikan perbaikan signifikan dengan mengestimasi parameter dispersi dan menghasilkan QAIC yang lebih kecil dibanding model binomial. Namun demikian, model ini tetap mempertahankan struktur varians linier terhadap \(\mu(1-\mu)\), yang secara empiris belum sepenuhnya mencerminkan pola variabilitas data leaf blotch.
Model quasi-custom menggunakan variance function berbentuk kuadrat,
\[ V(\mu) = (\mu(1-\mu))^2, \]
yang sesuai dengan formulasi Wedderburn dan interpretasi yang diberikan oleh Millar (2011). Hasil estimasi menunjukkan bahwa nilai dispersion mendekati satu (\(\approx 0.989\)) dan residual deviance hampir sama dengan derajat bebas model. Kondisi ini mengindikasikan bahwa struktur varians yang diasumsikan telah menangkap variasi data secara memadai.
Walaupun nilai QAIC quasi-binomial sedikit lebih kecil dibanding quasi-custom, pemilihan model dalam kerangka quasi-likelihood tidak hanya didasarkan pada kriteria informasi, tetapi juga pada kesesuaian struktur varians terhadap karakteristik data. Dalam konteks ini, model quasi-custom lebih konsisten secara teoretis dan empiris karena mampu menggambarkan peningkatan variabilitas yang lebih tajam terhadap mean.
Dengan mempertimbangkan:
maka model yang paling tepat untuk data leaf blotch adalah:
\[ \textbf{Quasi-Custom dengan } V(\mu)=(\mu(1-\mu))^2. \]
Model ini memberikan representasi varians yang lebih realistis, menghasilkan inferensi yang stabil, dan lebih sesuai dengan struktur empiris data dibandingkan model binomial maupun quasi-binomial standar.
G. Referensi
- Faraway, J. (2025). faraway: Functions and datasets for books by Julian Faraway (R package version 2.0.0) [PDF manual]. The Comprehensive R Archive Network (CRAN). https://cran.r-project.org/web/packages/faraway/faraway.pdf
- McCullagh, P., & Nelder, J. A. (1989). Generalized linear models (2nd ed.). Chapman & Hall
- Millar, R. B. (2011). Maximum likelihood estimation and inference: With examples in R, SAS and ADMB. John Wiley & Sons, Ltd.
- Wedderburn, R. W. M. (1974). Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method. Biometrika, 61(3), 439–447. https://doi.org/10.2307/2334725