backPain                 package:gnm                 R Documentation

_D_a_t_a _o_n _B_a_c_k _P_a_i_n _P_r_o_g_n_o_s_i_s, _f_r_o_m _A_n_d_e_r_s_o_n (_1_9_8_4)

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

     Data from a study of patients suffering from back pain. Prognostic
     variables were recorded at presentation and progress was
     categorised three weeks after treatment.

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

     data(backPain)

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

     A data frame with 101 observations on the following 4 variables.

     _x_1 length of previous attack.

     _x_2 pain change.

     _x_3 lordosis.

     _p_a_i_n an ordered factor describing the progress of each patient
          with levels 'worse' < 'same' < 'slight.improvement' <
          'moderate.improvement' < 'marked.improvement' <
          'complete.relief'. 

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

     <URL: http://ideas.repec.org/c/boc/bocode/s419001.html>

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

     Anderson, J. A. (1984) Regression and Ordered Categorical
     Variables. _J. R. Statist. Soc. B_, *46(1)*, 1-30.

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

     set.seed(1)
     data(backPain)

     ### Re-express as count data
     library(nnet)
     .incidence <- class.ind(backPain$pain)
     .counts <- as.vector(t(.incidence))
     .rowID <- factor(t(row(.incidence)))
     backPain <- backPain[.rowID, ]
     backPain$pain <- C(factor(rep(levels(backPain$pain), nrow(.incidence)),
                               levels = levels(backPain$pain), ordered = TRUE),
                        treatment)

     ### Fit models described in Table 5 of Anderson (1984)

     ### Logistic family models
     noRelationship <- gnm(.counts ~ pain, eliminate = ~ .rowID,
                           family = "poisson", data = backPain)

     ## stereotype model
     oneDimensional <- update(noRelationship,
                              ~ . + Mult(pain - 1, x1 + x2 + x3 - 1),
                              iterStart = 3)

     threeDimensional <- update(noRelationship, ~ . + pain:(x1 + x2 + x3))

     ### Models to determine distinguishability in stereotype model
     .pain <- backPain$pain

     levels(.pain)[2:3] <- paste(levels(.pain)[2:3], collapse = " | ")
     fiveGroups <- update(noRelationship,
                          ~ . + Mult(as.ordered(.pain) - 1,
                                     x1 + x2 + x3 - 1))

     levels(.pain)[4:5] <- paste(levels(.pain)[4:5], collapse = " | ")
     fourGroups <- update(fiveGroups)

     levels(.pain)[2:3] <- paste(levels(.pain)[2:3], collapse = " | ")
     threeGroups <- update(fourGroups)

     ### Grouped continuous model, aka proportional odds model
     library(MASS)
     sixCategories <- polr(pain ~ x1 + x2 + x3, data = backPain)

     ### Obtain number of parameters and log-likelihoods for equivalent
     ### multinomial models as presented in Anderson (1984)
     logLikMultinom <- function(model){
         object <- get(model)
         if (inherits(object, "gnm")) {
             l <- logLik(object) + object$eliminate
             c(nParameters = attr(l, "df") - object$eliminate, logLikelihood = l)
         }
         else
             c(nParameters = object$edf, logLikelihood = -deviance(object)/2)
     }
     models <- c("threeDimensional", "oneDimensional", "noRelationship",
                 "fiveGroups", "fourGroups", "threeGroups", "sixCategories")
     t(sapply(models, logLikMultinom))

