wedderburn                package:gnm                R Documentation

_W_e_d_d_e_r_b_u_r_n _Q_u_a_s_i-_l_i_k_e_l_i_h_o_o_d _F_a_m_i_l_y

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

     Creates a 'link{family}' object for use with 'glm', 'gnm', etc.,
     for the variance function  [mu(1-mu)]^2 introduced by Wedderburn
     (1974) for response values in [0,1].

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

     wedderburn(link = "logit")

_A_r_g_u_m_e_n_t_s:

    link: The name of a link function.  Allowed are "logit", "probit"
          and "cloglog". 

_V_a_l_u_e:

     An object of class 'family'.

_N_o_t_e:

     The reported deviance involves an arbitrary constant (see
     McCullagh and Nelder, 1989, p330); for estimating dispersion, use
     the Pearson chi-squared statistic instead.

_A_u_t_h_o_r(_s):

     David Firth

_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.

_S_e_e _A_l_s_o:

     'glm', 'gnm', 'family'

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

     set.seed(1)
     data(barley)  ##  data from Wedderburn (1974), see ?barley

     ###  Fit Wedderburn's logit model with variance proportional to the
     ###  square of mu(1-mu)
     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(LETTERS[1:9], as.character(1:9), "X"))

     ###  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)

