R/select.parfm.R

  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
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
################################################################################
#  Computation of AIC and BIC of many models of class 'parfm'                  #
################################################################################
#                                                                              #
#  The function 'select.parfm' computes the AIC and BIC values                 #
#    of parametric frailty models with                                         #
#    different baseline hazards and                                            #
#    different frailty distributions                                           #
#                                                                              #
#  Its parameters are                                                          #
#   - formula  : a formula object, with the response                           #
#                on the left of a ~ operator,                                  #
#                and the terms on the right.                                   #
#                The response must be a survival object                        #
#                as returned by the Surv function.                             #
#   - cluster  : the name of the variable in data containing cluster IDs       #
#   - strata   : the name of the variable in data containing strata IDs        #
#   - data     : a data.frame in which to interpret the variables named        #
#                in the formula.                                               #
#   - dist     : the vector of the names of the baseline hazards               #
#   - frailty  : the vector of the names of the frailty distribution           #
#   - method   : the optimization method (See optim())                         #
#   - maxit    : the maximum number of iterations (See optim())                #
#   - showtime : show the execution time of each model? (See parfm())          #
#   - correct  : the correction to use in case of many                         #
#                events per cluster to get finite likelihood values.           #
#                When correct!=0 the likelihood is divided by                  #
#                10 ^ (#clusters * correct) for computation,                   #
#                but the value of the log-likelihood in the output             #
#                is the re-adjusted value.                                     #
#                It is used only for models with Positive Stable frailty       #
#                                                                              #
#                                                                              #
#                                                                              #
#  The function returns a list with elements                                   #
#   - AIC : a table with AIC values of the required models                     #
#           with one line   per baseline hazard distribution and               #
#           with one column per frailty         distribution                   #
#   - BIC : a table with BIC values of the required models                     #
#           with one line   per baseline hazard distribution and               #
#           with one column per frailty         distribution                   #
#                                                                              #
#                                                                              #
#   Date: December 21, 2011                                                    #
#   Last modification on: December 2, 2016                                     #
################################################################################

select.parfm <- function(formula,
                         cluster=NULL,
                         strata=NULL,
                         data,
                         inip=NULL,
                         iniFpar=NULL,
                         dist=c("exponential",
                                "weibull",
                                "inweibull",
                                "frechet",
                                "gompertz",
                                "loglogistic",
                                "lognormal",
                                "logskewnormal"),
                         frailty=c("none",
                                   "gamma",
                                   "ingau",
                                   "possta",
                                   "lognormal"),
                         method="BFGS",
                         maxit=500,
                         Fparscale=1,
                         correct=0){
  warn <- getOption("warn")
  options(warn = -1)
  
  dist <- match.arg(dist, several.ok = TRUE)
  if (('inweibull' %in% dist) & 'frechet' %in% dist) {
      warning("Baselines 'inweibull' and 'frechet' are synonyms.")
      dist <- setdiff(dist, 'frechet')
  }
  
  res <- list(AIC = NULL, BIC = NULL)
  res$AIC <- res$BIC <- matrix(NA, length(dist), length(frailty),
                               dimnames = list(dist, substr(frailty, 1, 6)))
  cat(paste("\n\n### - Parametric frailty models - ###",
            "Progress status:",
            "  'ok' = converged",
            "  'nc' = not converged\n",
            "                Frailty",
            "Baseline           ",
            sep = "\n"))
  cat(c(none = " none  ", gamma = " gamma ", ingau = " invGau", 
        possta = " posSta", lognormal = " lognor")[frailty])
  for (d in dist) {
    cat("\n")
    cat(c(exponential = "exponential.......",
          weibull     = "Weibull...........",
          inweibull     = "inWeibull.........",
          frechet     = "inWeibull...........",
          gompertz    = "Gompertz..........",
          loglogistic = "loglogistic.......",
          lognormal   = "lognormal.........",
          logskewnormal   = "logskewnormal.....")[d])
    for (f in frailty) {
      cat("..")
      model <- try(parfm(formula = formula, 
                         cluster = cluster,
                         strata = strata,
                         data = data,
                         inip = inip,
                         iniFpar = iniFpar,
                         dist = d,
                         frailty = f,
                         method = method,
                         maxit = maxit,
                         Fparscale = Fparscale,
                         showtime = FALSE,
                         correct = correct),
                   silent = TRUE)
      if (!("try-error" %in% class(model))) {
        res$AIC[d, substr(f, 1, 6)] <- AIC(model)
        res$BIC[d, substr(f, 1, 6)] <- BIC(model)
        cat("ok....")
      } else cat("nc....")
    }
  }
  cat("\n")
  class(res) <- "select.parfm"
  
  options(warn = warn)
  return(res)
}

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.