1Icens <- function(first.well, last.well, first.ill, formula, 2 model.type=c("MRR","AER"), breaks, 3 boot=FALSE, alpha=0.05, keep.sample=FALSE, data) 4{ 5 ## Create follow-up matrix containing three event times 6 fu.expression <- substitute(cbind(first.well, last.well, first.ill)) 7 fu <- if (missing(data)) { 8 eval(fu.expression) 9 } 10 else { 11 eval(fu.expression, data) 12 } 13 ## Check consistency of arguments 14 missing.f1 <- is.na(fu[,1]) 15 missing.f2 <- is.na(fu[,2]) 16 if (any(missing.f1 & missing.f2)) { 17 stop("You must supply at least one of \"first.well\" and \"last.well\"") 18 } 19 if (any(fu[,1] > fu[,2], na.rm=TRUE) | any(fu[,2] > fu[,3], na.rm=TRUE)) { 20 stop("Some units do not meet: first.well < last.well < first.ill" ) 21 } 22 ## Fill in any gaps 23 fu[,1][missing.f1] <- fu[,2][missing.f1] 24 fu[,2][missing.f2] <- fu[,1][missing.f2] 25 ## Recensor cases that fall after the last break point 26 is.censored <- fu[,3] > max(breaks) 27 is.censored[is.na(is.censored)] <- FALSE 28 fu[is.censored,3] <- NA 29 30 exp.dat <- expand.data(fu, formula, breaks, data) 31 32 model.type <- match.arg(model.type) 33 if (missing(formula)) { 34 fit.out <- with(exp.dat, fit.baseline(y, rates.frame)) 35 lambda <- coef(fit.out) 36 } 37 else { 38 fit.out <- switch(model.type, 39 "MRR"=with(exp.dat, fit.mult(y, rates.frame, cov.frame)), 40 "AER"=with(exp.dat, fit.add(y, rates.frame, cov.frame))) 41 lambda <- coef(fit.out$rates) 42 } 43 44 beta <- if (is.null(fit.out$cov)) numeric(0) else coef(fit.out$cov) 45 params <- c(lambda,beta) 46 if (boot) { 47 nboot <- ifelse (is.numeric(boot), boot, 100) 48 boot.coef <- matrix(NA, nrow=nboot, ncol=length(lambda) + length(beta)) 49 colnames(boot.coef) <- names(params) 50 51 for (i in 1:nboot) { 52 subsample <- sample(nrow(fu), replace=TRUE) 53 exp.dat <- expand.data(fu[subsample,], formula, breaks, data[subsample,]) 54 if (missing(formula)) { 55 sim.out <- with(exp.dat, fit.baseline(y, rates.frame, params)) 56 boot.coef[i,] <- coef(sim.out) 57 } 58 else { 59 sim.out <- switch(model.type, 60 "MRR"=with(exp.dat, 61 fit.mult(y, rates.frame, cov.frame, params)), 62 "AER"=with(exp.dat, 63 fit.add(y, rates.frame, cov.frame, params))) 64 boot.coef[i,] <- switch(model.type, 65 "MRR"=c(coef(sim.out[[1]]), 66 coef(sim.out[[2]])), 67 "AER"=coef(sim.out[[1]])) 68 } 69 } 70 ci.quantiles=c(0.5, alpha/2, 1 - alpha/2) 71 boot.ci <- t(apply(boot.coef,2,quantile,ci.quantiles)) 72 73 lower.ci.lab <- 74 paste("lower ", formatC(100 * alpha/2, format="f", digits=1),"%", sep="") 75 upper.ci.lab <- 76 paste("upper ", formatC(100 * (1-alpha/2), format="f", digits=1),"%", 77 sep="") 78 79 colnames(boot.ci) <- c("median", lower.ci.lab, upper.ci.lab) 80 fit.out$boot.ci <- boot.ci 81 if (keep.sample) { 82 fit.out$sample <- boot.coef 83 } 84 } 85 86 class( fit.out ) <- "Icens" 87 attr( fit.out, "model" ) <- model.type 88 return( fit.out ) 89} 90