parfm: R/anova.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
################################################################################
#  anova method for class 'parfm'                                              #
################################################################################
#                                                                              #
#  anova.parfm implemets the 'anova' function the objects of class 'parfm'     #
#                                                                              #
#                                                                              #
#   Date: June 5, 2012                                                         #
#   Last modification on: June 5, 2012                                         #
################################################################################

anova.parfm <- function(object,
                        ...) {
  if (length(list(object, ...)) > 1)
    return(anova.parfmlist(list(object, ...)))
  else{
    if (!inherits(object, "parfm"))
      stop(paste("The argument must be a parametric frailty model,",
                 "object of class 'parfm'"))
    
    termlist <- attr(object, "terms")
    
    nmodels <- length(termlist)
    df <- integer(nmodels + 1)
    loglik <- double(nmodels + 1)
    df[1] <- 0
    mycall <- attributes(object)$call
    mycall[["showtime"]] <- FALSE
    mycall[[2]] <- update.formula(mycall[[2]], . ~ 1)
    loglik[1] <- attr(eval(mycall), "loglik")
    for (i in 1:nmodels) {
      df[1 + i] <- df[i] + 1
      mycall[[2]] <- eval(parse(text=paste(
        "update.formula(mycall[[2]], . ~ . +", termlist[i], ")")))
      loglik[1 + i] <- attr(eval(mycall), "loglik")
    }
    table <- data.frame(loglik=loglik, 
                        Chisq=c(NA, 2*diff(loglik)),
                        Df=c(NA, diff(df)))
    table[['Pr(>|Chi|)']] <- 1- pchisq(table$Chisq, table$Df)
    row.names(table) <- c('NULL', termlist)
    title <- paste("Analysis of Deviance Table\n",
                   "Parametric frailty model: response is ",
                   deparse(attr(object, "call")[
                     match("formula", names(attr(object, "call")))][[1]][[2]]),
                   "\nTerms added sequentially (first to last)\n", 
                   sep = "")
    structure(table, heading = title, class = c("anova", "data.frame"))    
  }
}

################################################################################
#  anova method for a list of objects of class 'parfm'                         #
################################################################################
#                                                                              #
#                                                                              #
#   Date: June 5, 2012                                                         #
#   Last modification on: June 5, 2012                                         #
################################################################################

anova.parfmlist <- function(object) {
  if (!is.list(object)) stop("First argument must be a list")
  if (!all(unlist(lapply(object, function(x) inherits(x, 'parfm')))))
    stop("Argument must be a list of 'parfm' models")
  
  nmodels <- length(object)
  if (nmodels == 1) # only one model remains
    return(anova.parfm(object[[1]]))
  
  ### *** Controls *** #########################################################
  responses <- as.character(unlist(lapply(object, function(x) {
    deparse(formula(attr(x, "call"))[[2]])})))
  sameresp <- (responses == responses[1])
  if (!all(sameresp)) {
    object <- object[sameresp]
    warning(paste("Models with response", deparse(responses[!sameresp]), 
                  "removed because response differs from", "model 1"))
  }
  
  whichargs <- names(attr(object[[1]], "call")[
    -match(c("", "formula"), names(attr(object[[1]], "call")))])
  for (arg in whichargs) {
    for (i in 2:nmodels) {
      if (attr(object[[i]], "call")[arg][[1]] != 
        attr(object[[1]], "call")[arg][[1]]) stop(paste(
          "All models must have the same value for argument '",
          arg, "'", sep=""))
    }
  }
  
  ##############################################################################
  
  loglik <- df <- double(nmodels)
  for (i in 1:nmodels) {
    loglik[i] <- attr(object[[i]], "loglik")
    df[i] <- sum(!is.na(object[[i]][, "ESTIMATE"]))
  }
  
  table <- data.frame(loglik, Chisq= c(NA, abs(2*diff(loglik))), 
                      Df= abs(c(NA, diff(df))))
  
  tfun <- function(x) paste(as.character(delete.response(
    terms(formula(attr(x, "call"))))), collapse=' ')
  
  variables <- lapply(object, tfun)
  dimnames(table) <- list(1:nmodels, 
                          c("loglik", "Chisq", "Df"))
  title <- paste("Analysis of Deviance Table\n Parametric frailty model: response is ",
                 responses[1]) 
  topnote <- paste(" Model ", format(1:nmodels), ": ", variables, 
                   sep = "", collapse = "\n")
  table[['P(>|Chi|)']] <- 1-pchisq(table$Chisq, table$Df)
  
  structure(table, heading = c(title, topnote), 
            class = c("anova", "data.frame"))
}

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.