polf                  package:VGAM                  R Documentation

_P_o_i_s_s_o_n-_O_r_d_i_n_a_l _L_i_n_k _F_u_n_c_t_i_o_n

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

     Computes the Poisson-ordinal transformation, including its inverse
     and the first two derivatives.

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

     polf(theta, earg = stop("'earg' must be given"), inverse = FALSE,
          deriv = 0, short = TRUE, tag = FALSE)

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

   theta: Numeric or character. See below for further details.

    earg: Extra argument for passing in additional information. This
          must be list with component 'cutpoint'. If 'polf()' is used
          as the link function in 'cumulative' then one should choose
          'reverse=TRUE, parallel=TRUE, intercept.apply=TRUE'.

 inverse: Logical. If 'TRUE' the inverse function is computed.

   deriv: Order of the derivative. Integer with value 0, 1 or 2.

   short: Used for labelling the 'blurb' slot of a 'vglmff-class'
          object.

     tag: Used for labelling the linear/additive predictor in the
          'initialize' slot of a 'vglmff-class' object. Contains a
          little more information if 'TRUE'.

_D_e_t_a_i_l_s:

     The Poisson-ordinal link function (POLF) can be applied to a 
     parameter lying in the unit interval. Its purpose is to link
     cumulative probabilities associated with an ordinal response
     coming from an underlying Poisson distribution. 

     The arguments 'short' and 'tag' are used only if 'theta' is
     character.

     See 'Links' for general information about 'VGAM' link functions.

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

     See Yee (2006) for details.

_W_a_r_n_i_n_g:

     Prediction may not work on 'vglm' or 'vgam' etc. objects if this
     link function is used.

_N_o_t_e:

     Numerical values of 'theta' too close to 0 or 1 or out of range
     result in large positive or negative values, or maybe 0 depending
     on the arguments. Although measures have been taken to handle
     cases where 'theta' is too close to 1 or 0, numerical
     instabilities may still arise.

     In terms of the threshold approach with cumulative probabilities
     for an ordinal response this link function corresponds to the
     Poisson distribution (see 'poissonff') that has been recorded as
     an ordinal response using known cutpoints.

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

     Thomas W. Yee

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

     Yee, T. W. (2006) _Link functions for ordinal count data_, 
     (submitted for publication).

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

     'Links', 'poissonff', 'nbolf', 'golf', 'cumulative'.

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

     polf("prob", short=FALSE)
     polf("prob", tag=TRUE)

     p = seq(0.01, 0.99, by=0.01)
     earg = list(cutpoint=2)
     y = polf(p, earg=earg)
     y. = polf(p, earg=earg, deriv=1)
     max(abs(polf(y, earg=earg, inv=TRUE) - p)) # Should be 0

     ## Not run: 
     par(mfrow=c(2,1), las=1)
     plot(p, y, type="l", col="blue", main="polf()")
     abline(h=0, v=0.5, col="red", lty="dashed")

     plot(p, y., type="l", col="blue",
          main="(Reciprocal of) first POLF derivative")
     ## End(Not run)

     # Rutherford and Geiger data
     ruge = c(57,203,383,525,532,408,273,139,45,27,10,4,0,1,1)
     y = 0:14
     yy = rep(y, times=ruge)
     length(yy)  # 2608 1/8-minute intervals
     cutpoint = 5
     yy01 = ifelse(yy <= cutpoint, 0, 1)
     earg = list(cutpoint=cutpoint)
     fit = vglm(yy01 ~ 1, binomialff(link="polf", earg=earg))
     coef(fit, matrix=TRUE)
     exp(coef(fit))

     # Another example
     nn = 1000
     x2 = sort(runif(nn))
     x3 = runif(nn)
     mymu = exp( 3 + 1 * x2 - 2 * x3)
     y1 = rpois(nn, lambda=mymu)
     cutpoints = c(-Inf, 10, 20, Inf)
     cuty = Cut(y1, breaks=cutpoints)
     ## Not run: 
     plot(x2, x3, col=cuty, pch=as.character(cuty))
     ## End(Not run)
     table(cuty) / sum(table(cuty))
     fit = vglm(cuty ~ x2 + x3, fam = cumulative(link="polf",
                reverse=TRUE, parallel=TRUE, intercept.apply=TRUE,
                mv=TRUE, earg=list(cutpoint=cutpoints[2:3])),
                trace=TRUE)
     fit@y[1:5,]
     fitted(fit)[1:5,]
     predict(fit)[1:5,]
     coef(fit)
     coef(fit, matrix=TRUE)
     constraints(fit)
     fit@misc$earg

