TABLE OF CONTENTS
methodsBiv/bivrec.formula [ Functions ]
NAME
bivrec.formula --- method for formulas
FUNCTION
This function is the one dispatched to by bivrec, when a formula is used as the response. Its function is to read the input, ensure the response is a Surv2 object, evaluate the model, and create an Anderson-Gill data frame in the format required by bivrec.agdata, then dispatch to bivrec.agdata for the fitting.
SYNOPSIS
92 bivrec.formula <- function(formula, data = parent.frame(), ...)
INPUTS
See package documentation
OUTPUTS
An object of class bivrec. See package documentation
SOURCE
95 { 96 # in part based on coxph function 97 # evaluate model in parent frame and get covariates, terms and response 98 call <- match.call() 99 m <- match.call(expand.dots = FALSE) 100 if(is.matrix(eval(m$data, sys.parent()))) m$data <- as.data.frame(data) 101 m$...<-NULL 102 m[[1]] <- as.name("model.frame") 103 special <- c("id", "cluster", "strata") 104 Terms <- if (missing(data)) terms(formula, special) else 105 terms(formula, special, data = data) 106 m$formula <- Terms 107 oldNAoption <- getOption("na.action"); options(na.action = na.fail) 108 m <- eval(m, sys.parent()) 109 options(na.action = oldNAoption) 110 n <- nrow(m) 111 112 # Extract response and check that it is in order 113 resp <- model.extract(m, "response") 114 if (!is.Surv2(resp)) stop("model response must be a Surv2 object") 115 start <- resp[, "start"] 116 stop <- resp[, "stop"] 117 delta <- resp[, "status1"] 118 Delta <- resp[, "status2"] 119 processnames <- attr(resp, "processnames") 120 dropx <- NULL 121 clusterind <- attr(Terms, "specials")$cluster 122 123 # Cluster handling 124 clusternames <- NULL 125 if(length(clusterind) > 0){ 126 cluster <- m[, clusterind] 127 if(is.factor(cluster)) clusternames <- levels(cluster) 128 if(is.numeric(cluster)) clusternames <- as.character(unique(cluster)) 129 i <- as.numeric(as.factor(cluster)) 130 tempc <- untangle.specials(Terms, "cluster", 1:10) 131 ord <- attr(Terms, "order")[tempc$terms] 132 if (any(ord > 1)) stop("Cluster cannot be used in an interaction") 133 dropx <- c(dropx, tempc$terms) 134 }else{ 135 cluster <- rep(1, n) 136 i <- rep(1, n) 137 } 138 139 # ID handling 140 idind <- attr(Terms, "specials")$id 141 if(length(idind) == 0) stop("a subject id (unique within clusters) is required") 142 else{ 143 id <- m[, idind] 144 subjnames <- paste(cluster, id, sep = ".") 145 j <- rep(0, n) 146 if(is.factor(id)) levels(id) <- 1:length(levels(id)) 147 for(ii in 1:max(i)){ 148 iids <- unique(id[i == ii]) 149 jj <- 1 150 for(thisid in iids){ 151 j[i == ii & id == thisid] <- jj 152 jj <- jj + 1 153 } 154 } 155 tempi <- untangle.specials(Terms, "id", 1:10) 156 ord <- attr(Terms, "order")[tempi$terms] 157 if (any(ord > 1)) stop("id cannot be used in an interaction") 158 dropx <- c(dropx, tempi$terms) 159 } 160 161 # Stratum handling 162 stratind <- attr(Terms, "specials")$strata 163 if(length(stratind) > 0) 164 { 165 strat <- m[, stratind] 166 ustrat <- unique(strat) 167 stratnames <- as.character(ustrat) 168 r <- rep(0, n) 169 for(rr in 1:length(ustrat)){ 170 r[strat == ustrat[rr]] <- rr 171 } 172 tempr <- untangle.specials(Terms, "strata", 1:10) 173 ord <- attr(Terms, "order")[tempr$terms] 174 if (any(ord > 1)) stop("strata cannot be used in an interaction") 175 dropx <- c(dropx, tempr$terms) 176 }else{ 177 r <- rep(1, n) 178 stratnames <- "1" 179 } 180 181 # Compute event counter 182 Ki <- table(i * 1e6 + j) 183 k <- unlist(as.vector(sapply(Ki, function(x) 1:x))) 184 185 # drop the specials and construct the model matrix 186 newTerms <- if(length(dropx)) Terms[ - dropx] else Terms 187 X <- model.matrix(newTerms, m) 188 X <- X[, -1, drop = FALSE] 189 190 # Construct the data frame in the format of bivrec.agdata 191 agdata <- as.data.frame(cbind(i, j, k, r, start, stop, delta, Delta, X)) 192 sortord <- order(agdata$i, agdata$j) 193 agdata <- agdata[sortord, ] 194 subjnames <- unique(subjnames[sortord]) 195 196 # Basic check that the input is in order 197 check <- checkinput.bivrec(agdata, clusternames, subjnames, 198 stratnames, processnames) 199 200 # Pass it to bivrec.agdata 201 fit <- bivrec(agdata, ...) 202 fit$call <- call 203 204 # Clean up the output 205 fit <- cleanbivrecoutput(fit, clusternames, subjnames, stratnames, processnames) 206 return(fit) 207 }