Module: Quasi likelihood

Paper: Regression Analysis III

Module: Quasi likelihood

Principal investigator: Dr. Bhaswati Ganguli,Professor, Department of Statistics, University of Calcutta

Paper co-ordinator: Dr. Bhaswati Ganguli,Professor, Department of Statistics, University of Calcutta

Content writer: Sayantee Jana, Graduate student, Department of Mathematics and Statistics, McMaster University Sujit Kumar Ray,Analytics professional, Kolkata

Content reviewer: Department of Statistics, University of Calcutta

I f(y) =e

yθ−b(θ) φ +c(y,φ)

, takinga(φ) =φ

I E(Y) =µ=b⁰(θ)

I V ar(Y) =b⁰⁰(θ)φ=V(µ)φ

I f(y) =e

yθ−b(θ) φ +c(y,φ)

, takinga(φ) =φ

I E(Y) =µ=b⁰(θ)

I V ar(Y) =b⁰⁰(θ)φ=V(µ)φ

I f(y) =e

yθ−b(θ) φ +c(y,φ)

, takinga(φ) =φ

I E(Y) =µ=b⁰(θ)

I V ar(Y) =b⁰⁰(θ)φ=V(µ)φ

I L= ^yθ−b(θ)_φ +c(y, φ)

I Differentiating Lw.r.tµ E[δL

δµ] = E[δL δθ.δθ

δµ]

= E[δL δθ]δθ

δµ

= E[y−b⁰(θ) φ ]δθ

δµ

I L= ^yθ−b(θ)_φ +c(y, φ)

I Differentiating Lw.r.tµ E[δL

δµ] = E[δL δθ.δθ

δµ]

= E[δL δθ]δθ

δµ

= E[y−b⁰(θ) φ ]δθ

δµ

I L= ^yθ−b(θ)_φ +c(y, φ)

I Differentiating Lw.r.tµ E[δL

δµ] = E[δL δθ.δθ

δµ]

= E[δL δθ]δθ

δµ

= E[y−b⁰(θ) φ ]δθ

δµ

I V ar[^δL_δµ] =V ar[^δL_δθ](_δµ^δθ)² =E[−^δ_δθ^2L₂](^δθ_δµ)²

I Where

E[δL

δθ]² = E[−δ²L δθ²]

= b⁰⁰(θ) φ

1 [b⁰⁰(θ)]²

= 1

b⁰⁰(θ).φ

I V ar[^δL_δµ] =V ar[^δL_δθ](_δµ^δθ)² =E[−^δ_δθ^2L₂](^δθ_δµ)²

I Where

E[δL

δθ]² = E[−δ²L δθ²]

= b⁰⁰(θ) φ

1 [b⁰⁰(θ)]²

= 1

b⁰⁰(θ).φ

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

I Very often the underlying distribution of the population may not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis. Recourse is then taken to a quasi - likelihood function which mimics the likelihood function and satisfied (I) and (II).

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

I Very often the underlying distribution of the population may not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis. Recourse is then taken to a quasi - likelihood function which mimics the likelihood function and satisfied (I) and (II).

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

I Very often the underlying distribution of the population may not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis. Recourse is then taken to a quasi - likelihood function which mimics the likelihood function and satisfied (I) and (II).

I Defineq = _V^y−µ_(µ).φ

I E(q) = 0and V(q) = _V(µ).φ¹

I q satisfied (I) and (II) and hence it can be used as a proxy to

δL δµ

I Defineq = _V^y−µ_(µ).φ

I E(q) = 0and V(q) = _V(µ).φ¹

I q satisfied (I) and (II) and hence it can be used as a proxy to

δL δµ

I Defineq = _V^y−µ_(µ).φ

I E(q) = 0and V(q) = _V(µ).φ¹

I q satisfied (I) and (II) and hence it can be used as a proxy to

δL δµ

I Defineq = _V^y−µ_(µ).φ

I E(q) = 0and V(q) = _V(µ).φ¹

I q satisfied (I) and (II) and hence it can be used as a proxy to

δL δµ

I Q=Rµ y

y−u φV(u)du

I δQ

δµ = _φV^y−µ_(µ) =q

I Q as a substitute forL

I Only requirement to write q or Qis to identifyV(µ) as a function of µ.

I To estimate to β

¯ : maximizePn

i=1Qi w.r.tβ

¯ .

I Q=Rµ y

y−u φV(u)du

I δQ

δµ = _φV^y−µ_(µ) =q

I Q as a substitute forL

I Only requirement to write q or Qis to identifyV(µ) as a function of µ.

I To estimate to β

¯ : maximizePn

i=1Qi w.r.tβ

¯ .

I Q=Rµ y

y−u φV(u)du

I δQ

δµ = _φV^y−µ_(µ) =q

I Q as a substitute forL

I Only requirement to write q or Qis to identifyV(µ) as a function of µ.

I To estimate to β

¯ : maximizePn

i=1Qi w.r.tβ

¯ .

I Q=Rµ y

y−u φV(u)du

I δQ

δµ = _φV^y−µ_(µ) =q

I Q as a substitute forL

I Only requirement to write q or Qis to identifyV(µ) as a function of µ.

I To estimate to β

¯ : maximizePn

i=1Qi w.r.tβ

¯ .

I Q=Rµ y

y−u φV(u)du

I δQ

δµ = _φV^y−µ_(µ) =q

I Q as a substitute forL

I Only requirement to write q or Qis to identifyV(µ) as a function of µ.

I To estimate to β

¯ : maximizePn

i=1Qi w.r.tβ

¯ .

u^∗_r = δ δβ_r

n

X

i=1

Q_i

=

n

X

i=1

δQ_i δµi

δµ_i δηi

δη_i δβr

=

n

X

i=1

q_i(δµi

δη_i)x_ir

=

n

X

i=1

y_i−µ_i V(µi)φ(δµ_i

δηi

)xir

Q = Z µ

y

y−u φu du

= 1

φ[ylogu−u]^µ_y

= 1

φ[ylogµ−µ−ylogy+y]

u^∗ = 1Xⁿ yi−µi

(δµi

)x

Q = Z µ

y

y−u φu du

= 1

φ[ylogu−u]^µ_y

= 1

φ[ylogµ−µ−ylogy+y]

u^∗ = 1Xⁿ yi−µi

(δµi

)x

Q = Z µ

y

y−u φu du

= 1

φ[ylogu−u]^µ_y

= 1

φ[ylogµ−µ−ylogy+y]

u^∗ = 1Xⁿ yi−µi

(δµi

)x

Q = Z _µ

y

y−u φ du

= 1

φ[u.y−u² 2 ]^µ_y

= 1

φ[µy−µ²

2 −y²+y² 2 ]

= 1

2φ[−y²+ 2µy−µ²]

= −(y−µ)²

I V(µ) =µ(1−µ)

Q = Z µ

y

y−u φ(u(1−u))du

= 1 φ{

Z µ y

y

u(1−u)du− Z µ

y

du 1−u}

= 1 φ{y

Z µ y

1−u+u u(1−u)du−

Z µ y

du 1−u}

= 1 φ{y

Z µ y

1 udu+y

Z µ y

du 1−u −

Z µ y

du 1−u}

= 1

φ{[ylogu]^µ_y −[ylog(1−u)]^µ_y + [log(1−u)]^µ_y}

= 1

φ{ylogµ−ylogy−ylog(1−µ)+ylog(1−y)+log(1−µ)−log(1−y)}

= 1

φ{ylog µ

1−µ+log(1−µ)−ylog y

1−y −log(1−y)}

ylog_1−µ^µ +log(1−µ) ←log likelihood for the Binomial distribution

Q=

n

X

i=1

Q_i Where

Qi = Z µi

yi

yi−u φV(u)du D= 2

n

X

i=1

Z yi

µi

y_i−u V(u) du

I estimation : quasi score vector : u

¯= ^δQ

δβ

n ¯

=X

i=1

y_i−µ_i φV(µ_i)

δµ_i δβ

¯

= 1

φDV⁻¹(y

¯

−µ

¯) Where D= ((^δµi

δβ

¯ )) V =diag((V(µi)))

¯y−µ

¯ =











Y1−µ1

Y₂−µ₂ ...

−











A = E[− δ²Q δβ

¯δβ

¯

0]

= −1 φ

n

X

i=1

E[(y_i−µ_i) δ δβ

¯

0

1 V(µ_i)(δµ_i

δβ

¯

) + 1 V(µ_i)(δµ_i

δβ

¯ ) δ

δβ

¯

0(y_i−µ_i)]

= 1 φ

n

X

i=1

(δµ_i δβ

¯ ) 1

V(µi)(δµ_i δβ

¯

0)]

= 1

φD⁰V⁻¹D

At them^th iteration βˆ

¯

(m) = βˆ

¯

(m−1)

+ [A⁻¹u

¯] β

¯

=βˆ

¯

(m−1)

= βˆ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯

−µ

¯)]_{β= ˆˆ β}(m−1)

At them^th iteration βˆ

¯

(m) = βˆ

¯

(m−1)

+ [A⁻¹u

¯] β

¯

=βˆ

¯

(m−1)

= βˆ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯

−µ

¯)]_{β= ˆˆ β}(m−1)

I E(ˆβ

¯) =β

¯+O(n⁻¹) D(ˆβ

¯) =φ(D⁰V⁻¹D)⁻¹+O(n⁻¹)

φ(D⁰V⁻¹D)⁻¹ ← involves the overdispersion parameters

I Estimate of φ: No estimate of φcan be directly obtained from the Quasi - likelihood.

Conventionally φis estimated as - φˆ= 1

n−p

n

X

i=1

(y_i−µˆ_i)² V(ˆµ_i) This is moment estimation.

I E(ˆβ

¯) =β

¯+O(n⁻¹) D(ˆβ

¯) =φ(D⁰V⁻¹D)⁻¹+O(n⁻¹)

φ(D⁰V⁻¹D)⁻¹ ← involves the overdispersion parameters

I Estimate of φ: No estimate of φcan be directly obtained from the Quasi - likelihood.

Conventionally φis estimated as - φˆ= 1

n−p

n

X

i=1

(y_i−µˆ_i)² V(ˆµ_i) This is moment estimation.

not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis.

I To overcome this a quasi - likelihood function which mimics the likelihood function and satisfies (I) and (II) is used.

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis.

I To overcome this a quasi - likelihood function which mimics the likelihood function and satisfies (I) and (II) is used.

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis.

I To overcome this a quasi - likelihood function which mimics the likelihood function and satisfies (I) and (II) is used.

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis.

I To overcome this a quasi - likelihood function which mimics the likelihood function and satisfies (I) and (II) is used.

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

not be known (e.g in case of overdispersed data). In such cases , the likelihood can not be obtained and hence inferences can not be drawn on its basis.

I To overcome this a quasi - likelihood function which mimics the likelihood function and satisfies (I) and (II) is used.

I E[^δL_δµ]= 0 ...(I)

I V ar[^δL_δµ]= ¹

b⁰⁰(θ).φ ...(II)

I Quasi deviance function : D= 2Pn i=1

Ryi

µi

yi−u V(u)du

I estimation : quasi score vector : u

¯= ¹_φDV⁻¹(y

¯

−µ

¯)

I Quasi information matrix : A= _φ¹D⁰V⁻¹D

I Estimator : βˆ

¯

(m)= ˆβ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯−µ

¯)]_{β= ˆˆ β}(m−1)

I Estimate of φ: φˆ= _n−p¹ Pn i=1

(yi−ˆµi)² V(ˆµi)

I Quasi deviance function : D= 2Pn i=1

Ryi

µi

yi−u V(u)du

I estimation : quasi score vector : u

¯= ¹_φDV⁻¹(y

¯

−µ

¯)

I Quasi information matrix : A= _φ¹D⁰V⁻¹D

I Estimator : βˆ

¯

(m)= ˆβ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯−µ

¯)]_{β= ˆˆ β}(m−1)

I Estimate of φ: φˆ= _n−p¹ Pn i=1

(yi−ˆµi)² V(ˆµi)

I Quasi deviance function : D= 2Pn i=1

Ryi

µi

yi−u V(u)du

I estimation : quasi score vector : u

¯= ¹_φDV⁻¹(y

¯

−µ

¯)

I Quasi information matrix : A= _φ¹D⁰V⁻¹D

I Estimator : βˆ

¯

(m)= ˆβ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯−µ

¯)]_{β= ˆˆ β}(m−1)

I Estimate of φ: φˆ= _n−p¹ Pn i=1

(yi−ˆµi)² V(ˆµi)

I Quasi deviance function : D= 2Pn i=1

Ryi

µi

yi−u V(u)du

I estimation : quasi score vector : u

¯= ¹_φDV⁻¹(y

¯

−µ

¯)

I Quasi information matrix : A= _φ¹D⁰V⁻¹D

I Estimator : βˆ

¯

(m)= ˆβ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯−µ

¯)]_{β= ˆˆ β}(m−1)

I Estimate of φ: φˆ= _n−p¹ Pn i=1

(yi−ˆµi)² V(ˆµi)

I Quasi deviance function : D= 2Pn i=1

Ryi

µi

yi−u V(u)du

I estimation : quasi score vector : u

¯= ¹_φDV⁻¹(y

¯

−µ

¯)

I Quasi information matrix : A= _φ¹D⁰V⁻¹D

I Estimator : βˆ

¯

(m)= ˆβ

¯

(m−1)

+ [(D⁰V⁻¹D)⁻¹D⁰V⁻¹(y

¯−µ

¯)]_{β= ˆˆ β}(m−1)

I Estimate of φ: φˆ= _n−p¹ Pn i=1

(yi−ˆµi)² V(ˆµi)

# install.packages("qcc")

require(qcc) ## contains the function qcc.overdispersion.test

# data from Wetherill and Brown (1991) pp. 212--213, 216--218:

x <- c(12,11,18,11,10,16,9,11,14,15,11,9,10,13,12, 8,12,13,10,12,13,16,12,18,16,10,16,10,12,14) size <- rep(50, length(x))

qcc.overdispersion.test(x,size)

x <- c(11,8,13,11,13,17,25,23,11,16,9,15,10,16,12, 8,9,15,4,12,12,12,15,17,14,17,12,12,7,16) qcc.overdispersion.test(x)

data(breslow.dat, package = "robust") names(breslow.dat)

head(breslow.dat)

summary(breslow.dat[c(6, 7, 8, 10)])

# plot distribution of post-treatment seizure counts opar <- par(no.readonly = TRUE)

par(mfrow = c(1, 2)) attach(breslow.dat)

hist(sumY, breaks = 20, xlab = "Seizure Count", main = "Distribution of Seizures")

boxplot(sumY ~ Trt, xlab="Treatment",main="Group Comparisons") par(opar)

family = poisson()) summary(fit)

# interpreting model parameters coef(fit)

exp(coef(fit))

# evaluating overdispersion

qcc.overdispersion.test(breslow.dat$sumY, type = "poisson")

# fit model with quasipoisson

fit.od <- glm(sumY ~ Base + Age + Trt, data = breslow.dat, family = quasipoisson())