wedderburn.barley: Percent of barley leaves affected by leaf blotch

Description Format Details Source References Examples

Description

Percent of leaf area affected by leaf blotch on 10 varieties of barley at 9 sites.

Format

A data frame with 90 observations on the following 3 variables.

y

Percent of leaf area affected, 0-100.

site

Site factor, 9 levels

gen

Variety factor, 10 levels

Details

Incidence of Rhynchosporium secalis (leaf blotch) on the leaves of 10 varieties of barley grown at 9 sites in 1965.

Source

Wedderburn, R W M (1974). Quasilikelihood functions, generalized linear models and the Gauss-Newton method. Biometrika, 61, 439–47. http://www.jstor.org/stable/2334725.

Wedderburn credits the original data to an unpublished thesis by J. F. Jenkyn.

References

McCullagh, P and Nelder, J A (1989). Generalized Linear Models (2nd ed).

R. B. Millar. Maximum Likelihood Estimation and Inference: With Examples in R, SAS and ADMB. Chapter 8.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
data(wedderburn.barley)
dat <- wedderburn.barley
dat$y <- dat$y/100

require(lattice)
dotplot(gen~y|site, dat, main="wedderburn.barley")

# Use the variance function mu(1-mu).  McCullagh page 330
# Note, 'binomial' gives same results as 'quasibinomial', but also a warning
m1 <- glm(y ~ gen + site, data=dat, family="quasibinomial")
summary(m1)

# Same shape (different scale) as McCullagh fig 9.1a
plot(m1, which=1, main="wedderburn.barley")

# Compare data and model
dat$pbin <- predict(m1, type="response")
dotplot(gen~pbin+y|site, dat, main="wedderburn.barley: observed/predicted")


# Wedderburn suggested variance function: mu^2 * (1-mu)^2
# Millar shows how to do this explicitly.
wedder <- 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.99,y)) }),
             name="(mu(1-mu))^2")
m2 <- glm(y ~ gen + site, data=dat, family=quasi(link="logit", variance=wedder))

## Not run: 
# Alternatively, the 'gnm' package has the 'wedderburn' family.
require(gnm)
m3 <- glm(y ~ gen + site, data=dat, family="wedderburn")
summary(m3)
# Similar to McCullagh fig 9.2
plot(m3, which=1)

# Compare data and model
dat$pwed <- predict(m3, type="response")
dotplot(gen~pwed+y|site, dat)

## End(Not run)

agridat documentation built on May 30, 2017, 5:11 a.m.