Contoh_quasi_likelihood
Dariani Matualage
2025-02-25
R Code yang digunakan diambil dari buku (Millar, 2011, hal 178-180)
library(writexl)
library(readxl)
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(tidyr)
library(gnm)Input data
Data yang digunakan adalah data proporsi dauh yang terserang bercak (Millar, 2011 hal 177). Data diukur pada 10 varitas Barley di setiap 9 lokasi.
BarleyTable=matrix(barley$y,ncol=10,
dimnames=list(paste("site",1:9,sep=""),paste("vrty",1:10,sep="")));BarleyTable## vrty1 vrty2 vrty3 vrty4 vrty5 vrty6 vrty7 vrty8 vrty9 vrty10
## site1 0.0005 0.0150 0.0100 0.2000 0.250 0.080 0.050 0.05 0.050 0.250
## site2 0.0000 0.0000 0.1270 0.3750 0.550 0.165 0.050 0.50 0.050 0.425
## site3 0.0000 0.0005 0.0125 0.2625 0.050 0.295 0.100 0.10 0.250 0.500
## site4 0.0010 0.0005 0.0125 0.0250 0.400 0.200 0.050 0.50 0.750 0.375
## site5 0.0025 0.0030 0.0250 0.0050 0.055 0.435 0.500 0.25 0.500 0.950
## site6 0.0005 0.0075 0.1660 0.0001 0.010 0.010 0.750 0.50 0.750 0.625
## site7 0.0050 0.0030 0.0250 0.0300 0.060 0.050 0.050 0.75 0.750 0.950
## site8 0.0130 0.0300 0.0250 0.0250 0.011 0.050 0.001 0.05 0.750 0.950
## site9 0.0150 0.0750 0.0000 0.0001 0.025 0.050 0.050 0.10 0.175 0.950Quasi-Binomial
Pendugaan dilakukan dengan menggunakan model Quasi-likelihood
quasibin.fit=glm(y~variety+site,family=quasibinomial,data=barley)Hasil pendugaan
summary(quasibin.fit)##
## Call:
## glm(formula = y ~ variety + site, family = quasibinomial, data = barley)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.0546 1.4219 -5.665 2.84e-07 ***
## variety2 0.1501 0.7237 0.207 0.836289
## variety3 0.6895 0.6724 1.025 0.308587
## variety4 1.0482 0.6494 1.614 0.110910
## variety5 1.6147 0.6257 2.581 0.011895 *
## variety6 2.3712 0.6090 3.893 0.000219 ***
## variety7 2.5705 0.6065 4.238 6.58e-05 ***
## variety8 3.3420 0.6015 5.556 4.39e-07 ***
## variety9 3.5000 0.6013 5.820 1.51e-07 ***
## varietyX 4.2530 0.6042 7.039 9.38e-10 ***
## siteB 1.6391 1.4433 1.136 0.259870
## siteC 3.3265 1.3492 2.466 0.016066 *
## siteD 3.5822 1.3444 2.664 0.009510 **
## siteE 3.5831 1.3444 2.665 0.009493 **
## siteF 3.8933 1.3402 2.905 0.004875 **
## siteG 4.7300 1.3348 3.544 0.000697 ***
## siteH 5.5227 1.3346 4.138 9.38e-05 ***
## siteI 6.7946 1.3407 5.068 3.00e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.08877777)
##
## Null deviance: 40.803 on 89 degrees of freedom
## Residual deviance: 6.126 on 72 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 8Koefisien dugaan dari model quasi-binomial tersebut identik dengan koefisien yang diperoleh dari regresi logistik binomial biasa. namun, penduga parameter dispersi yang diperoleh adalah 0.0888. Artinya bahwa model quasi-binomial mengasumsikan var (Yij) = 0.0888 pij(1-pij) dan galat baku telah diskalakan dengan faktor (0.0888)^(1/2) dibandingkan dengan yang diperoleh dari regresi logistik binomial.
membuat plot deviance residuals vs linear predictors
#Use which=1 for a plot of deviance resids vs linear predictors
plot(quasibin.fit,which=1)
plot residual dari model quasi-binomial menunjukkan peningkatan yang jelas pada besarnya sisaan dengan peningkatan nilai prediktor liniernya, diikuti dengan kemungkinan penurunan untuk nilai prediktor linier yang lebih dari 0. ini menunjukkan bahwa mungkin ada lebih banyak keragaman dalam data untuk nilai prediktor linier mendekati 0 (yang sesuai dengan p mendekati 0.5) daripada untuk nilai ekstrim prediktor linier (yang sesuai dengan p mendekati 0 atau 1).
Menggunakan Wedderburn (1974):
BarleyVar = 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" )model pendugaan
wedderburn.fit=glm(y~site+variety,
family=quasi(link="logit",variance=BarleyVar),
data=barley)
plot(wedderburn.fit,which=1) # deviance residuals vs linear predictors
Hasil pendugaan
summary(wedderburn.fit)##
## Call:
## glm(formula = y ~ site + variety, family = quasi(link = "logit",
## variance = BarleyVar), data = barley)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.92238 0.44465 -17.817 < 2e-16 ***
## siteB 1.38312 0.44465 3.111 0.00268 **
## siteC 3.86006 0.44465 8.681 8.20e-13 ***
## siteD 3.55700 0.44465 8.000 1.54e-11 ***
## siteE 4.10786 0.44465 9.239 7.53e-14 ***
## siteF 4.30536 0.44465 9.683 1.13e-14 ***
## siteG 4.91810 0.44465 11.061 < 2e-16 ***
## siteH 5.69489 0.44465 12.808 < 2e-16 ***
## siteI 7.06763 0.44465 15.895 < 2e-16 ***
## variety2 -0.46735 0.46870 -0.997 0.32204
## variety3 0.07881 0.46870 0.168 0.86695
## variety4 0.95408 0.46870 2.036 0.04547 *
## variety5 1.35263 0.46870 2.886 0.00515 **
## variety6 1.32854 0.46870 2.835 0.00595 **
## variety7 2.34007 0.46870 4.993 4.01e-06 ***
## variety8 3.26258 0.46870 6.961 1.30e-09 ***
## variety9 3.13549 0.46870 6.690 4.10e-09 ***
## varietyX 3.88727 0.46870 8.294 4.34e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasi family taken to be 0.9885464)
##
## Null deviance: 252.155 on 89 degrees of freedom
## Residual deviance: 71.175 on 72 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 20plot sisaan dari model quasi-likelihood Wedderburn menunjukkan peningkatan yang besar. deviasi sisaan sebesar 71.175 adalah statistik chi-square Person yang telah didefinisikan dalam BarleyVar, yang secara implisit mengasumsikan prediktor linier=1. deviasi sisaan ini secara kebetulan memberikan nilai yang sangat dekat dengan derajat bebas 72.
galat baku yang dihasilkan mengasumsikan parameter dispersi sebesar 71.175/72 ~ 0.989.