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