barley                  package:gnm                  R Documentation

_J_e_n_k_y_n'_s _D_a_t_a _o_n _L_e_a_f-_b_l_o_t_c_h _o_n _B_a_r_l_e_y

_D_e_s_c_r_i_p_t_i_o_n:

     Incidence of _R. secalis_ on the leaves of ten varieties of barley
     grown at nine sites.

_U_s_a_g_e:

     data(barley)

_F_o_r_m_a_t:

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

     _y the proportion of leaf affected (values in [0,1])

     _s_i_t_e a factor with 9 levels 'A' to 'I'

     _v_a_r_i_e_t_y a factor with 10 levels 'c(1:9, "X")'

_N_o_t_e:

     This dataset was used in Wedderburn's original paper (1974) on 
     quasi-likelihood.

_S_o_u_r_c_e:

     Originally in an unpublished Aberystwyth PhD thesis by J F Jenkyn.

_R_e_f_e_r_e_n_c_e_s:

     Gabriel, K R (1998).  Generalised bilinear regression. 
     _Biometrika_  *85*, 689-700.

     McCullagh, P and Nelder, J A (1989)  _Generalized Linear Models_
     (2nd ed).   Chapman and Hall.

     Wedderburn, R W M (1974).  Quasilikelihood functions, generalized
     linear models  and the Gauss-Newton method.  _Biometrika_ *61*,
     439-47.

_E_x_a_m_p_l_e_s:

     data(barley)
     set.seed(1)

     ##  Fit Wedderburn's logit model with variance proportional to [mu(1-mu)]^2
     logitModel <- glm(y ~ site + variety, family = wedderburn, data = barley)
     fit <- fitted(logitModel)
     print(sum((barley$y - fit)^2 / (fit * (1-fit))^2))
     ##  Agrees with the chi-squared value reported in McCullagh and Nelder 
     ##  (1989, p331), which differs slightly from Wedderburn's reported value.

     ##  Fit the biplot model as in Gabriel (1998, p694)
     biplotModel <- gnm(y ~ -1 + Mult(site, variety, multiplicity = 2),
                        family = wedderburn, data = barley)
     barleySVD <- svd(matrix(biplotModel$predictors, 10, 9))
     A <- sweep(barleySVD$v, 2, sqrt(barleySVD$d), "*")[, 1:2]
     B <- sweep(barleySVD$u, 2, sqrt(barleySVD$d), "*")[, 1:2]
     ##  These are essentially A and B as in Gabriel (1998, p694), from which
     ##  the biplot is made by
     plot(rbind(A, B), pch = c(levels(barley$site), levels(barley$variety)))

     ##  Fit the double-additive model as in Gabriel (1998, p697)
     variety.binary <- factor(match(barley$variety, c(2,3,6), nomatch = 0) > 0,
                             labels = c("rest", "2,3,6"))
     doubleAdditive <- gnm(y ~ variety + Mult(site, variety.binary),
                           family = wedderburn, data = barley)
     ##  It is unclear why Gabriel's chi-squared statistics differ slightly
     ##  from the ones produced in these fits.  Possibly Gabriel adjusted the
     ##  data somehow prior to fitting?

